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1 . INTRODUCTION 


The objective of this contract is to develop a computer pro- 
gram to predict the inviscid, transonic flow field about isolated 
nacelles. Furthermore, the problem is to be formulated to solve 
Euler's equations without any approximation (such as small dis- 
turbances) and hence the terminology "exact" solution. 

The flow field is complicated by the presence of imbedded 
shock waves, an engine-inlet interface and exhaust plumes. Fur- 
thermore, the transonic nacelles of interest have a very slender 
but blunt cowl lip. This creates two distinct length scales, the 
length of the nacelle and the cowl lip radius that can differ by 
several orders of magnitude. These aspects of the flow field 
present many numerical difficulties. 

Our approach to the problem is to calculate the nacelle flow 
field using the method of time -dependent computations (TDC) . Al- 
though at the time of the issuance of this contract, other approaches 
to transonic flow calculations existed, we felt that TDC offered 
the most effective means of meeting the goals of the contract. 

This report is divided into two major sections: Section 2, 

the Final Report and Section 3, the Computer Program Documentation. 
The Final Report contains a discussion of our approach, the mathe- 
matical formulation of the problem, some results and conclusions. 

The Computer Program Documentation contains one section geared 
toward the engineering user of the program and a second section 
for programming considerations . 
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2. FINAL REPORT 

2.1 GENERAL FEATURES 

2.1.1 Approach 

Our basic approach to solving the inviscid transonic flow 
field about an isolated nacelle consists of using the method of 
time-dependent computations (TDC) . A time-dependent flow can be 
analyzed by solving an initial and boundary value problem, that 
is well suited for numerical solution. Steady state solutions 
are obtained from asymptotic, large time results. This technique 
has been used successfully in the past in a number of relatively 
simple transonic problems such as the supersonic blunt body flow 
(Ref. 1) , the flow in a choked nozzle (Ref. 2), and the flow past 
a boattail (Ref. 3). 

Another feature of our approach is the handling of imbedded 
shock waves. We feel that the most accurate and efficient pro- 
cedure is to consider all shocks as discontinuities, across which 
the Rankine-Hugoniot relations must be satisfied. Successful com- 
putations of this type have been made for one dimensional flows 
(Ref. 4) and for a few selected two dimensional transonic problems 
(Ref. 3) . 

The details of our basic approach are discussed in Section 2.2. 
Many of these techniques have been evolved from the work described 
in Refs. 3 and 5. Other specific features have been developed by 
necessity throughout the course of this effort, which will be dis- 
cussed in the next section. 

2.1.2 History 

Our first attempt at solving the transonic nacelle problem 
consisted in using the approach in Refs. 6 and 7 . Namely, a blunt 
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nacelle was considered to accelerate from rest to the free stream 
Mach number in a specified time interval. The computation would 
be performed between the body surface and the perturbation wave 
front, shown schematically in Fig. 2.1. Several problems appeared 
using this formulation. Firstly, as the wave front moved further 
from the body, the mesh resolution deteriorated severely. It was 
necessary to maintain sufficient resolution immediately behind the 
wave front which caused a further loss in resolution near the body 
surface. Secondly, we had a polar toroidal coordinate system 
adjacent to a cylindrical coordinate system. The interface between 
these two systems had to be treated very carefully to avoid error. 
These problems were greatly amplified when we began to treat 
nacelles approaching the leading edge radii of the test case, 
r/L 0.05-0.010. For these cases there was no way to maintain 
resolution within the limitations of the storage of the computer. 
Furthermore, it was observed that the perturbation front was ex- 
tremely weak and, thus, not essential for the computation. 

Therefore, since it would be necessary to have the program 
handle thin-nosed nacelles, the first approach had to be dropped. 
While problems were developing with the original computation pro- 
cedure, an attempt was made to compute nacelles as if they had an 
infinitely thin, cusped-cowl lip. This approach, shown schemati- 
cally in Fig. 2.2 considered the entire infinite domain surrounding 
the nacelle. The basic methodology used here was an extension of 
Ref. 3. The results for the subsonic, flow-through nacelle were 
encouraging . 

The next step was to modify the cusped nose nacelle with an 
actual blunt nose. Here it was decided to use the same coordinate 
system stretchings of the cusped nose case with a small overlapped 
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Fig. 2.1 Coordinate System — First Approach 
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Fig. 2.2 Cusped - Nacelle Formulation 
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region of very fine mesh point resolution near the nose as in 
Fig. 2.3. The nose computation was carried on independently of 
the computation external to the nose with spacial and temporal 
interpolation at the interface boundaries. This procedure al- 
lowed us to have very fine resolution near the nose with a cor- 
responding small time step and a more sparse mesh outside the nose 
region with a correspondingly large time step. This approach, 
that was successfully implemented for the very slender nosed 
nacelles used in the experiments, was adopted as our "final" ap- 
proach. It is described in detail in Section 2.2. 

The final approach has been applied successfully to several 
test cases with a subsonic free stream and the results are de- 
scribed in Section 2.3. When higher free stream Mach numbers or 
lower inlet mass flow ratios were attempted, problems developed.. 
Namely, during the transient phase of the computation the flow 
goes supersonic inside the cowl lip. This situation appears to 
be valid on physical grounds. A shock should form inside the cowl 
lip and move upstream eventually "popping" out of the inlet. The 
stagnation streamline would then move inside to cowl lip and the 
flow would be subsonic inside the inlet. However, this process 
causes severe numerical problems. The shock must move upstream 
through the inlet past the cowl lip and away from the nacelle. It 
must pass through the overlapped nose region. The complexities of 
nose computation seem to preclude the representation of this 
"starting" shock as a discontinuity. 

Although we generally feel that it is desirable to represent 
shocks as discontinuities on the basis of accuracy and efficiency 
factors, we did make some attempts at "smearing" this shock. As 
a first step, we tried to change the entire formulation of the 
equations to conservation form. A preliminary study of using this 
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Fig. 2.3 Schematic Mesh Point Distribution 
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approach was performed on a simplified isolated boattail and the 
results were encouraging. However, when the procedure was used 
in the complete nacelle, computation problems developed. The 
calculation deteriorated in the vicinity of the nose region in- 
terface. This situation existed even without the presence of 
any imbedded shock waves. Apparently, the nose region interpola- 
tion procedure did not work properly with the conservation form 
variables so the approach was dropped. 

In another attempt to handle this starting shock, we con- 
sidered a crude smearing procedure. The points where the flow 
decelerated from supersonic to subsonic Mach numbers were deter- 
mined. The last supersonic point was computed using upwind dif- 
ferences. Since the shock could not pop out of the inlet, the 
appropriate value of the mass flow was never achieved. An ex- 
ample of this case is discussed in Section 2.3. 

Even if problems with the starting shock could be solved, 
there would be other difficulties in the transonic regime. Shocks 
also appear on the cowl outer surface. During the transient, at 
low mass flow ratios, these shocks may start right at the cowl lip 
and progress around the nose and downstream until the steady state 
position is obtained. Thus, similar difficulties to those for the 
starting shock are obtained. We have attempted to overcome the 
problems of the shock appearing at the nose by using the smearing 
technique just discussed. This procedure works adequately to delay 
the difficulties until the shock is located downstream of the nose. 
At this point, we feel the shock should be fit as a discontinuity. 
Once again, this situation presents a problem. Most of our ex- 
perience with imbedded shocks in two space dimensions (Ref. 3) 
consists of having the shock as a fixed internal boundary and 
having the mesh points move with the shock. However, complications 
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at the nose region seem to preclude this approach. The alterna- 
tive of having a shock move between node points has not been 
fully investigated and it appears to require further study before 
it can be applied to this problem. 

Our conclusions regarding the treatment of imbedded shocks 
for the transonic nacelle problem will be discussed in Section 2.4. 

2.2 MATHEMATICAL FORMULATION 

The basic aspects of the flow field for the nacelle are shown 
schematically in Fig. 2.4. Important features of the problem con- 
sist of: a long, thin nacelle with a blunt cowl lip of small 

radius of curvature, a blunt centerbody located inside the cowl, 
an inlet-engine interface where the mass flow is specified, an 
under -expanded supersonic jet plume and flow field boundaries ex- 
tending to infinity. The details of the formulation of these 
problems now follow. 

2.2.1 Basic Nacelle Computation 

The formulation of the basic nacelle computation includes the 
stretching of the flow field from the body surface to infinity and 
handling of the blunt cowl lip with a radius of curvature which is 
very small compared with the length of the nacelle. We feel that 
an efficient means of modeling the flow field is to first formulate 
the problem as if the nose of the nacelle were a cusp and then to 
modify the computation to include the effects of the blunt nose. 
Note that we are not approximating the nose shape, nor are we 
actually computing it as a^usp, but we are merely formulating the 
problem with a cusp and then modifying it. This feature will be- 
come clear when the nose region computation is discussed in Sec- 
tion 2.2.2. 
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Fig. 2.4 Transonic Nacelle — Flowfield Description 
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Consider the flow field for a cusp-nosed nacelle shown in 
Fig. 2.2. For convenience in the mapping functions that will be 
discussed later, the flow field is divided into seven regions. 
Regions 1 and 2 are upstream of the nacelle, Regions 3 and 4 are 
inside the inlet, Region 5 is the plume, Region 6 is downstream 
of the nacelle, and Region 7 is exterior to the nacelle. The 
origin of a cylindrical coordinate system, x,y, is located at 
the tip of the cusped cowl lip at the centerline. We wish to 
have a mesh point distribution in the physical plane that con- 
centrates mesh points near the nacelle surface and near the 
centerbody or axis. We also require the mesh points to be con- 
centrated near the leading and trailing edges of the nacelle, 
near the leading edge of the centerbody, and extending towards 
infinity (both upstream and downstream) with a sparse mesh point 
distribution. The appropriate mesh point distributions and 
stretchings to infinity can be accomplished using techniques dis- 
cussed in Refs. 3 and 6. We map each region to a computational 
square where uniform mesh spacing in the computational plane re- 
sults in the desired mesh point spacing in the physical plane. 

The stretching functions for each region are tabulated in Table 2-1. 
The parameters in the transformations x q , x-^, X£, x^, y Q , y 2 > 

(* 2 , may be selected by specifying the mesh spacing in the 

physical planes and will be explained in Section 3.1. 

The governing equations in the physical plane for the time- 
dependent inviscid flow over the nacelle may be written as 

P. + uP + vP + y(u + v + j ~) = 0 
t x y 'xy J y 


+ uu 

+ vu 

+ 

tP 

+ X (t) = 0 

X 

y 


X 

o 

+ uv 

+ w 

+ 

tP 

= 0 


t x y y 


( 1 ) 


S + uS + vS = 0 
t x y 
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G(Z;a 1 ,a 2 ) - ^ + a 2 log(Z) j log(Z) 








where P = log p, R = log p, S = P - 7R, t = exp^^^—^P + ^ sj , 
and where p,p are the pressure and density nondimensionalized by 
the free stream values, p ,p ; u and v are the velocities in 

00 00 

the x and y directions, respectively, nondimensionalized by 

V p /p and i=0 for two dimensional flow, and i = 1 for axi- 

00 00 

symmetric flow. The acceleration function X Q (t) will be dis- 
cussed at the end of this subsection in connection with initial 
conditions. 


The governing equations in the computational plane which are 
defined in Table 2-1 as 

X = X(x) 

Y = Y(x,y) (2) 

T = t 

now become 

p T + FP X + gp y + T(Au x + Bu y + Dv y + H) = 0 
u T + Fu x + Gu y + T(AP X + BP y ) + X Q = 0 

(3) 

v T + Fv x + Gv y + TDP y =0 
S T + FS X + GS y = 0 

where A=X , B = Y , D = Y , F=uA, G=uB+vD, H= j v/y, 
x x y 

and the derivatives of the stretching functions X x , Y x » Y^ are 
tabulated in Table 2-2. 

To solve Eqs . (3) numerically, a second order accurate finite 
difference approximation is used. Following the discussion of 
Ref. 8, a two level scheme used elsewhere by MacCormack (Ref. 9) 
is adopted. This two level scheme adapts to the regional makeup 
of our mesh point distribution. Namely, the flow field at points 
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Derivatives of Coordinate Stretchings 
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I 

along the computational interfaces labeled A, B, E, F, G in 
Fig. 2.2 may be evaluated in one region for the predictor stage 
and in the neighboring region for the corrector. For example, 
consider interface F. The predictor stage uses backward dif- 
ferences for the X-derivatives and the flow is determined along 
the interface from the equations in Region 1. The corrected 
values are then obtained from forward differences in Region 7. 
There will be no loss in accuracy provided the mesh spacings in 
the physical plane on each side of the interval are identical. 


A stable time step. At, is obtained by satisfying the 
Courant-Fredrichs-Lewy criterion. Namely, 

min (Ax, Ay) 

Ju 2 + V + V"t' 


At 


where 


A dX AV 

Ax = ax Ax 


Ay = || AX + §f AY 


The boundary conditions at the surface of the nacelle are 
the vanishing of the velocity normal to the wall. A procedure 
identical to that discussed in Ref. 3 is employed here. At the 
axis, the term H = j v/y is replaced by H = jOv^. 

Initially, at T = 0, the nacelle is assumed to be station- 
ary in a gas at rest; it accelerates to a transonic velocity in a 
finite period of time . The acceleration function is 



0 < T < 1/ co 
T > 1/co 


( 4 ) 
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where X Q is the abscissa of a reference point on the body with 
respect to a fixed frame. A more detailed discussion of the ini- 
tial conditions appears in Ref. 3. 

2.2.2 Nose Modifications 

Our technique for accurately computing the flow field near 
the blunt cowl consists of patching a small region over the cusped 
nacelle nose as shown in Fig. 2.5. The outer boundary of the nose 
region is a circle of radius, r = c . 

The computation procedure is as follows: 

1. Compute the entire flow field for the nacelle at 

t = t , , + At as if the nose were cusped for one 
new old 

time step. This time step. At, is the largest al- 
lowable, satisfying the C-F-L conditions, not in- 
cluding the nose region. 

2. Determine the largest allowable time step in the nose 

region. This time step will be smaller than the one 

used in step 1. Choose K -= At/At to be an 

r t nose 

integer. Typically, this value may be two or three. 

3. Determine the flow field at the outer boundary of the 
nose region by interpolation. 

4. Compute the flow field inside the nose region at 

t = t o ^ + A^Qggj using a polar coordinate system 
centered about the nacelle nose . 

5. Repeat steps 3 and 4 until t = t . 

r r nose new 

6. Determinate the flow field at the mesh points of the 
cusped nose formulation that is interior to the nose 
region. 

7. Repeat steps 1-6 for each time step, At. 
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CUSPED 

NACELLE 

MESH 




The governing equations in the nose region are the basic 
Euler equations written in a polar-toroidal (r,6) coordinate 
system (see Fig. 2.5). A computational plane (X,Y) (unit square) 
is obtained with the following transformation 

v _ 6 - tt/2 


TT 

= rj; bjex 

Y C - b(e) 

T = t 


(5) 


where r = b(0) is the surface of the cowl lip and r = c (a 
circle) is the outer boundary of the nose region. The governing 
equations may now be written as 

Bu„ Dv„ 


P T + EP y + FP x +y 


. . u . X . * . . 

Au„H I-— +j 

Y r r r J 

2 


u sin 

JL 

+ V 

cos Q 

H + 

r 

sin 

8 


- 0 


u T + Eu y + Fu x - ~ + TAP y - X q cos 6*0 


11X7 rr •• 

v T + Ev y +, Fv x + 7~'+“ (BP X + DP y ) + X Q cos e = 0 


S T + ES y + FS X = 0 


( 6 ) 


where B * 1 /tt, A = l/(c-b), D = (Y - l)b # /(c -b), E = uA + (v/r)D, 
F = Bv/r , and 

7 0 2-D 

1 Axisymmetric 


and X q is the acceleration of the body that has been discussed 
in connection with the initial conditions. Here, u and v are 
the nondimensional velocities in the radial and azimuthal direc- 
tions, respectively. 
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The C-F-L conditions for the nose region may be written as 


AT 


min(rA8 , r) 


nose 


Vu2 


+ + *Ty t 


where 

rAP = irrAX 

Ar = (c - b)AY - (Y - l)b # 7rAX . 

We will now discuss the interpolations necessary for the nose 
region computation. First, the values of the flow field at the 
outer boundary of the nose region must be determined from the 
values in the cusped nacelle formulation. Three typical situations 
appear as shown in Fig. 2.6. Cases 1 and 2 correspond to mesh 
situations near r = c and case 3 corresponds to the boundary 
6 = -rr/2 or 6 = 3tt/ 2 . Cusp nacelle mesh points 1 and 2 are found 
as the nearest set of nodes at least one of which is located ex- 
terior to the nose region. We designate and ^ n ew 

values of any function at these node points at time t ^ and 


t = t . , + At. The values of f at t 
new old nose 


t ,, + nAt 
old nose 


are interpolated so that 


fi - f old + < n/ V< f n 


- f l , ) 
new old 


where n = 1, 2, . . . , K t and i = 1 or 2 corresponds to points 1 

and 2. Point a corresponds to the boundary value whose value at 

t = t ,, + nAt remains to be determined. Point c corresponds 

old nose c 

to the nearest mesh point interior to the outer boundary. The 

values of the flow field at point c can be computed since they 

only depend upon the flow at point a at a time t = t ^ + 

(n - l)At . Point b corresponds to a point exterior to the 
nose 

nose region between points 1 and 2 along 9=6 =6 in cases 1 

Si c 
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CASE 1 


CASE 2 


CASE 3 


Fig. 2.6 Nose Region Interface Interpolation 



Fig. 2.7 Nose Region Interior Interpolation 



and 2, and along r = r = r in case 3. For cases 1 and 2, the 

d c 

value of f at point a, at t = t , , + nAt ^ is 

r ’ old nose 

f a - f c + - f c > < 7 > 

DC 

where for case 1, 

/ ” y l s 

f _ f + (-2 . f ) 

b 1 ^2 ‘ y l' 2 1 

and for case 2, 

f = f T + (— — )(f 9 - f t ) . 

b 1 \X2 - x^/ 2 1 

Similarly, for case 3, 

£ a - f c + (r^)( f b - V < 8 > 

1 c 

where 

f b - f i ;< f 2 ■ f i» • 

Now all the nose region may be computed K steps until a time 

t = t . , + KAt is reached. The values of the cusped 

new old nose 

nacelle points interior to the nose region may now be updated by 
interpolation of the nose region points. The interpolation is 
illustrated in Fig. 2.7. Points 1, 2, 3, and 4 are the node 
points in the nose region surrounding the cusped nacelle point a. 
Points b and c are points defined at 6 - 6 ^ at r = = r^ 

and r = r^ = r^, respectively. By linear interpolation we ob- 
tain 
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f b f 2 + Ve 3 - 0 j^) ^ f 4 ‘ f 2^ 


f c " f i + (e 3 - 0 1 )^ f 3 “ f P 


and 


f - f . + h 


r - r_. r - r.. 
^ r l r 2 + ^3 * r 4 


< f c ' V 


(9) 


The entire cycle is repeated for successive time steps 


2.2.3 Engine “-Inlet Interface 


The mass flow going through the engine at the engine-inlet 
interface is a boundary condition on the nacelle computation. 

From this specified mass flow the values of the flow variables 
P, u, v, S along interface C of Fig. 2.2 must be found. One 
technique consists of using unsteady characteristic compatibility 
relations coupled with conservation of mass along streamlines. 

This technique turns out to be quite cumbersome in practice. In 
order to avoid this difficulty, we have devised the computational 
device diagrammed in Fig. 2.8. Here we assume a choked convergent- 
divergent nozzle follows the interface. This is not the actual ex- 
haust nozzle, but is merely a computational technique to specify 
the required mass flow through the interface by fixing the area of 

•k 

the throat. Namely, the throat radius r is found to be 




r 

Iy + i 


l + 




zii ± n 

20v - i) 


.j/2 


( 10 ) 


where m r is the mass flow ratio, H is the height of the cowl 
lip above the centerline, and j = 0 for 2-D and j =1 for 
axisymmetric flow. 
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Fig. 2.8 Computational Device to Determine Engine-Inlet Conditions 



Fig. 2.9 Exhaust Plume for an Underexpanded Nozzle 
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Once the stream is accelerated to a supersonic velocity, the 
computation may be terminated by extrapolating the flow in the 
supersonic region. This device allows us to determine all the 
fluid properties at the interface by only specifying the mass 
flow. The detail of the flow through this nozzle, (the computa- 
tional device) is not pertinent to the nacelle calculation. 

2.2.4 Plume Computation 

The major objective of including the plume calculation is to 
determine the effect of the plume shape on the nacelle surface 
pressure distribution. We feel that this goal can be most ef- 
fectively achieved using a quasi-steady approach. Namely, the 
supersonic under -expanded jet (the case of interest) is computed 
using steady state equations with the plume boundary pressure 
distribution specified and the shape of the plume boundary re- 
sulting from the computation. This new jet boundary is then used 
in the time -dependent nacelle calculations for a specified number 
of time steps (typically about 100). Then the resulting pressure 
distribution from the nacelle calculation is input to the plume 
program and a new plume boundary results. This procedure is re- 
peated until the plume boundary pressure distribution and shape 
reach steady values . This approach should converge since rela- 
tively large changes in the pressure distribution will produce 
only very minor changes in the plume shape. 

The basic features of an under -expanded supersonic exhaust 
plume are illustrated in Fig. 2.9. This computation considers 
the flow field downstream to near the Mach disk. At this point, 
the shape of the plume boundary is extrapolated. The intercept- 
ing shock is handled as a discontinuity satisfying the Rankine- 
Hugoniot relations. 
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2.2. 4.1 Equations, Transformation 

The governing inviscid, steady Euler equations are written 


as 


uP + wP + y(u + w + j — ) = 0 
r z x z z J r 


uu + wu + tP = 0 
r z r 


uw + ww + tP =0 
r z z 


uS + wS =0 
r z 


(ID 


where P = log(p/p ref ), R = log(p/p ref ), S = P - yR. 


JP + — sl and u and w are the velocities in 

Y •' Y J 

the r and z directions, respectively, nondimensionalized by 


t = exp 




V 


r/p. 


ref' 'ref 

Note here that a change in notation has been made for the 
plume calculation discussion. The computer program for the plume 
was written independently of the nacelle program and, hence, used 
somewhat different variable names. For consistency, we will now 
adopt the notation of the plume program. The following chart lists 
the major changes. 


Plume 


Nacelle 


z 

— ¥ 

X 

r 

- 

y 

u 

— * 

V 

w 


u 


Furthermore, the nondimensionalization of the variables is dif- 
ferent. For the nacelle, p r = p , p c = p , arid Q r = 

’ *ref oo ’ r ref r oo’ x ref 

Vp / p . For the plume, p r = p . _ , p n = p. and Q r = 

r oo k qq r ’ r ref jet r ref jet x ref 

V Pj e t / Pj e t » where the subscript jet corresponds to the conditions 
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at the nozzle exhaust plane. The relationship of Pj et >Pj et to 
p ,p and the chamber conditions are 

00 r oo 



where is the stagnation pressure ratio and is the stagna- 

tion temperature ratio of the jet to the free stream. 

In general, the plume flow field will contain an intercepting 
shock as shown in Fig. 2.9. Here the flow field is split into two 
regions, I and II. Region I is bounded by the axis of the jet b(I) 
and the imbedded shock c(I) . Region II is bounded by the imbedded 
shock b(II) « c(I) and the plume boundary c(II) . When there is 
no intercepting shock, there will only be one region bounded by 
b(I) , the plume axis, and c(I), the plume boundary. A computa- 
tional frame for either region is obtained with the following 
transformation 



(13) 

Z ** z 

The equations are then transformed and rearranged to give the 
Z-derivatives explicitly as 
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where 


P„ = 


w„ = 


u r 


1 r 2 
2 


w 

- YT 

wR^ 

- tR 

2 


w 

- YT 

1 

U„(l 


X v 

1 

C3 — — * 

w 

s (u: 

X 


( 14 ) 


u x (uX r + «X z ) - tP x X ] 


R-, 


. Yu 


P x (uX r + *X z ) + Y( Ux X r 4- W x X z + j 7 


R 2 w X^ uX r + wX z> + tP X X z 


The equations are differenced using the MacCormack scheme described 

in Section 2.2.1. A inarching procedure in the Z-direction can be 

2 

used, provided w > y-v or the Mach number in the z-direction 
is supersonic. 


2. 2. 4. 2 Characteristics 

The evaluation of the jet boundary and the imbedded shock re- 
quire the use of characteristics. In general, the characteristic 
direction and compatibility relations in the computational plane 
can be developed using standard methods. The characteristic direc- 
tion is 


A = X 


X z uw ± a^g 

X + 2 2 

r w - a 


where 


a = y t 


P 


= Vm 2 - 


o 2 2 

M 2 = 


and the compatibility equation is 


ua u ± w3 


i (3 t(P z + ^P x ) + « <S Z + M x ) J- -y— 


(15) 


w - a 

where 6 = u/w. The C-F-L rule for this case can be represented 


as 


AZ = min 


(¥) • 


2. 2. 4. 3 Jet Boundary 

The pertinent characteristic reaching the plume boundary from 
inside the jet is A + as shown in the sketch in Fig. 2.10. This 
boundary is a streamline so that c = u/w and X = 1, whereby 


X = X 


r\ 


w/ 


so that 


A + = X 


- e + 


2 2 
8w + pa 

2 2 
w - a 


The compatibility relation becomes 


e z - - * +e x + T 

W 


+ A+ p x) . H|_ - 

w - a 


(16) 


(17) 
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with P^ prescribed from the nacelle computation. The above 
equation is integrated to give the position of the plume boundary 


at 


new 


Z + dZ. 


2. 2. 4. 4 Imbedded Shock 

In general, the intercepting shock does not immediately begin 
at the nozzle lip. Its presence is detected at a position down- 
stream of the nozzle exit by monitoring the pressure distribution 
in the r-direction at each z location. When the maximum pres- 
sure gradient grows in three consecutive z steps, a shock is fit. 
Initially, the shock is taken to be coincident with the A charac- 
teristic at the point of maximum gradient. Namely, 


uw - a p 

C Z ■ - * ' T 2 

w - a 


(18) 


As the shock develops, the values on the low pressure side are 
computed from the two characteristics reaching the shock (see 


Fig. 2.11). The values at the two points * 
tained by interpolation. 


and 


are ob- 


where 


P* = P Q + e + (P 1 - P~) 


6 2 + 6 ( S l ** 6 2 ^ 


C N ~ r 2 ~ *2 &Z 
Ar + Az(A“ - A") 


(19) 


The two characteristic compatibility relations may be simultaneously 
solved to get P^ and 6^ in terms of P*> e * and the other flow 
quantities at the * points. 
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where 




= (tP) + 

2 + 
= (*0 


2 

ua / u + wg N 

r \ 2 2 / 

w - a 



E x = - (t£) 
E 2 = - (w 2 ) 

+ D i P * + V* 




( ..P , . r ^ 

V 2 2/ 

w - a 


Az + E^P* + E 2 e“ 


( 20 ) 


The flow properties on the high pressure side at the shock are 
found from one characteristic relation combined iteratively with 
the Rankine-Hugoniot relations. The pertinent characteristic 
reaching the shock from the high pressure side is A as shown in 
Fig. 2.12. The values at the * point are found by interpolating 
between points 1 and 2, whereby 

p * = P i + e ( P 2 " P P 
e* = e 1 + e(e 2 - ep 

where 

c„ - c - A.Az 
N 1 

6 = ■" . 

Ar + Az(A 2 - A^) 
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The resulting compatibility relation may be integrated between 
the shock and the * point to yield 


:n. 


Hi 


A 

P + _s . hi -i." - "P) Ar 

* Sr Pr (» 2 -a 2 ) 




( 21 ) 


HI 


where Hi corresponds to the high pressure side of the shock. 
The Rankine-Hugoniot relations may be written as 


where 


- (7 - 1) 


P = P 
N„ . N. 
Hi Lo 


+ log 


7 + 1 


~ 2 
u. 


2 “1 
M = — 

n i a i 


( 22 ) 


7T. 


N Lo ’ 


u^ is the normal component of the velocity to the shock on the 
low pressure side of the shock, and Lo refers to the low pres- 
sure side of the shock. 


U, 


V N 1 + V N 2 


Lo 


Lo 


N, = - 


vT 


N, 


+ c ' 


Vl + c 2 

z 


XT is also a function of c . We take u„ , v„ to 
n H1 z 2 2 


In addition, 6 

be the normal and tangential components, respectively, on the high 
pressure side of the shock. From the Rankine-Hugoniot relations 
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and 

whereby 


The above Eqs . (21), (22), and (23), are solved in a trial and 
error procedure to determine c z at the new shock point. 

An example of these procedures for a plume computation will 
be presented in Section 2.3. 

2.3 DISCUSSION OF RESULTS 

The applications of the numerical procedures presented in 
Section 2.2 are now discussed. The examples include three com- 
putations of axisymmetric nacelle flow fields about a given geome- 
try with different mass flow ratios and Mach numbers and one com- 
putation of an isolated boattail with a plume. 

The nacelle geometry for all cases considered was the forebody 
of the NACA 1-85-100 nacelle designated Inlet No. 8 of NASA LRC 
test 264 as sketched in Fig. 2.13. Although the actual nacelle 
was 54 in. long, only the first 18 in. were used in the calcu- 
lations with a straight pipe assumed to be following aft of the 
cowl. The computations did not include the centerbody as shown 
in Fig. 2.13 and included a "throat” within the inlet to specify 


u. 


(Y - 1)M 2 + 2 

n. 


u. 


(Y + 1)M' 


n. 


v 2 = V 1 


U N t N 2 “ W N t N 1 
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Fig. 2.13 Nacelle Geometry (N ACA 1-85-100) 
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the appropriate mass flow as discussed in Section 2.2.3. The 
geometry was approximated with a series of 13 cubic fits with 
values of radius and slope at intersection points agreeing with 
the experimental geometry. The geometry routine is discussed in 
Section 3. 

All the nacelle calculations used approximately 2150 node 
points in the entire field with approximately 80 points on the 
surface of the nacelle and 19 of these on the cowl lip. The de- 
tails of the mesh point distributions are found in the test case 
in Section 3 . 

The first case computed is the flow field about the nacelle 
geometry discussed above at M = 0.4 and m^ = 0.847. The re- 
sults of the calculation after 2000 time steps, T = 2.43, are 
shown in Fig. 2.14. Values of Cp versus x/L axial distance 
near the cowl lip are shown compared to experimental data from 
run 20 Test 264 NASA LRC . 

As a second example, the M = 0.7, m r = 0.8715 flow over the 

same nacelle geometry as the first case was computed. Figure 2.15 

shows the calculated values of C versus x near the cowl lip 

P 

after 1600 time steps corresponding to a nondimensional time, 

T = 2.94. The experimental data of point 47, run 20 of the Langley 
test are in good agreement with the numerical results. The details 
of the computation in the vicinity of the cowl lip are shown in the 
isobar plot, Fig. 2.16. 

Next we attempted to compute the M = 0.9, m r = 0.885 flow 
over the same nacelle. In this calculation, a shock formed inside 
the cowl lip. As discussed previously in Section 2.1, we attempted 
to "smear" this shock. However, the shock remained inside the inlet 
and did not "pop" out. The effective result was that the mass flow 
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EXPERIMENT: 
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Fig. 2.14 Nacelle Surface Pressure Distribution Cowl Lip 
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Fig. 2.15 Nacelle Surface Pressure Distribution Cowl Lip 
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going into the inlet was not sufficiently low. Although the mass 
flow specified by the engine-inlet interface was 0.885, the 
actual calculated mass flow entering the inlet was 0.930. This 
fact is in evidence in the computed surface pressure distribution 
shown in Fig. 2.17. Here, the shape of the C versus x curve 
agrees with the data but is displaced. This phenomenon corresponds 
to an incorrectly matched mass flow. 

As an example of the boattail plume computation, we considered 
the flow over a straight semi-infinite pipe at M^ = 0.7. The jet 
had an initial Mach number of 3.0 with a ratio of jet static 
pressure to free stream static pressure of 3.0. The results of 
the computed plume shape, imbedded shock location and Mach number 
distributions within the jet are shown in Fig. 2.18. The calcula- 
tion ran 500 steps to a time T = 12.4 with the plume shape re- 
vised every 50 steps. The number of plume iterations was more 
than sufficient with the plume shape changing less than two per- 
cent during the final iteration. 

2.4 CONCLUSIONS 

The basic objectives of this task have been to develop a com- 
puter program for predicting the inviscid transonic flow field 
about a nacelle. Furthermore, it was desired to obtain an ex- 
tremely accurate solution to this problem. Our previous experi- 
ence in calculating transonic flows over simple bodies dictated 
our approach of using the method of time -dependent computations 
(TDC) . However, as we have described in Sections 2.1, 2.2, and 
2.3, the inherent complexities of this problem coupled with our 
basic approach have prevented us from obtaining a satisfactory 
solution at this time. 
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Fig. 2.17 Nacelle Surface Pressure Distribution Cowl Lip 
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We do not feel that there is anything fundamentally wrong 
with our methods. Rather, the sum total of complex features due 
to geometry, engine-inlet interfaces, exhaust plume and coordinate 
stretchings have combined to form an extremely complicated program. 
This type of program cannot accommodate our shock fitting procedures 
and, hence, the difficulties at transonic Mach numbers. Nontheless, 
the numerical techniques developed here will find some applications 
for the calculation of subsonic nacelle flow fields. 

In order to solve the transonic nacelle problem at this time, 
we feel that another approach may be found more suitable. Relaxa- 
tion methods have proven quite effective in solving two and three 
dimensional transonic airfoil problems (e.g., Jameson's work) and 
more recently for axisymmetric boattails (e.g., J. South). How- 
ever, some important technical details have to be developed before 
a relaxation approach can be used for the nacelle problem. Some 
examples are the determination of a suitable mapping so that the 
nacelle surface becomes a coordinate surface at the transformed 
plane, and the development of procedures to specify the inlet mass 
flow and to handle the exhaust plume computation. 

At present we feel that TDC has limited applicability for 
complicated transonic flow calculations. Ultimately, h clever, for 
transonic flows with strong shock waves and for time-varying flow 
fields, time-dependent methods should prove to be quite useful . 
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3. COMPUTER PROGRAM DOCUMENTATION 

3.1 USER ORIENTED DOCUMENTATION 

In this section we discuss all the details necessary to . 
enable the engineering user to run the computer program for the 
time-dependent nacelle calculation. Firstly, the general fea- 
tures and over-all logic flow are presented. The program usage 
and operation are then explained in terms of a complete discus- 
sion of program inputs. Finally, the accuracy and limitations 
of the program will be described. 

3.1.1 General Features and Over-all Logic Flow 

Here we present a qualitative description of all features of 
the computer program pertinent to the engineering user. A dis- 
cussion of the options, nacelle geometry, mesh point distribution, 
plume geometry, output options, plume computation, and over-all 
logic flow now follow. 

3 . 1 . 1 . 1 Options 

The basic computer program for the time-dependent computation 
of the inviscid flow field about a nacelle is designated program 
15C . This program has several major options. Namely, it can 
handle either a complete nacelle or an isolated boattail. There 
are also options for including an exhaust plume and for specify- 
ing either two dimensional or axisymmetric flow. In addition, 
the program may be started from time T = 0 directly or from a 
tape input generated in a previous run. The implementation of 
these features will be discussed in the section dealing with in- 
put and output. Other critical factors necessary to run the pro- 
gram consist of the geometry and mesh point distribution. 



3. 1.1. 2 Nacelle Geometry 


The geometry routine supplied as part of the computer program, 
called SUBROUTINE WALL, is relatively simple to implement. The 
nacelle surface is divided into three areas, the cowl lip, ex- 
ternal cowl, and internal cowl, as shown in Fig. 3.1. Each of 
these areas is divided into an arbitrary number of sections. The 
user inputs the number of divisions for each region and the value 
of y and dy/dx at each division point. In addition, a parame- 
ter is input at each point denoting whether a cubic or a straight 
line is to be fit between successive divisions. The program auto- 
matically performs the appropriate curve fits. The details of the 
user geometry input now follow. 

Firstly, the length L of the nacelle, the radius of the cowl 

lip r n , and the height above the centerline of the center of the 

cowl lip, H. Note that all barred qualities discussed here are 

dimensional. All lengths in the program will then be nondimens ional 

ized with respect to L. The radius of the cowl lip is taken from 

the nacelle blueprint. It is not necessary for this lip to be 

circular and the value of r^ is somewhat arbitrary. However, the 

smaller the r^, the greater the resolution near the nose, that 

also corresponds to increased computer running time. Typically, 

for the NACA 1-185-100 inlet, we used an r = 0.2 in. as is shown 

n 

in Section 3.2.4. 

Next, the cowl lip must be considered. This region is divided 
into JNOS intervals as shown in Fig. 3.2. The values of x,y and 
dy/dx are input at each division. These values must be in order 
starting at J = 1 corresponding to x = r^ on the upper surface 
to J = JNOS which is x = r^ on the lower surface. The program 
automatically curve fits r as a cubic in 8 between each interval 
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Fig. 3.1 Nacelle Geometry Notation 



Fig. 3.2 Cowl Lip Notation 
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Note that Cartesian values of x, y, dy/dx are input, with the 
program automatically converting to polar coordinates. It is 
recommended that one interval should correspond to x = 0 (the 
value of dy/dx at this point is infinite and any finite input 
value will be accepted) . Generally three to five intervals in 
the nose region should be adequate. 

The next part of the geometry input is the external cowl. Here, 
there are JOUT divisional points as shown in Fig. 3.3. The first in- 
tersection J => 1 corresponds to the cusped nacelle point and is 
automatically input as x ■ - r Q + r n and y = H, dy/dx =0. A cubic 
fit is used between J «= 1 and J = 2 . From J = 2 to J = JOUT 
values of x, y, dy/dx must be input along with the value of the 
parameter LOUT(J) . LOUT(J) = .3 for a cubic fit between J and J + 1 
and LOUT ( J ) » 1 for a straight line between J and J + 1. Also it 
is necessary to input the actual number of intersection points which are 
to be input called JOUTB. Since J - 1 is specified in the program, 

JOUTB ■= J OUT-1 . Note that we define r «= 3r . 

on 

The last region to be considered is the internal cowl as de- 
picted in Fig. 3.4. Here we have JINB divisional points with 
x, y, and dy/dx input for each point. Also LIN(J) is input, 
which determines the type of curve fit to be used between J and 

J + 1. The values at J = 1 are preset with x = - r + r , 

on 

y = H, dy/dx = 0, and LIN(l) = 3. Furthermore, values at the 
next to last intersection point JIN-1 are predetermined from 
the mass flow specification (Section 2.2.3) and the last value of 
JIN where x = L, y = H, and dy/dx = 0. Thus, only values be- 
tween J = 2 and J = JIN-2 need to be input . The total number 
of input intersectional points JINB = JIN-3. 

When the boattail option of the program is used, only the data 
for the external cowl need to be input. The plume geometry input 
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J = JOUT 




will be discussed separately. The nacelle geometry routine does 
not explicitly have a centerbody, although provisions for this 
situation have been maintained throughout the computer program. 

The nacelle geometry inputs will be summarized in Section 3.1.2. 

3 . 1 . 1 . 3 Mesh Point Distribution and Coordinate Stretchings 

Regions 1, 2, 3, 4, 6, and 7 shown in Fig. 2.2 have NC(LREG) 
mesh points in the x-direction and MC(LREG) mesh points in the 
y-direction. Region 8 has MC(8) mesh points in the r-direction 
and NC(8) in the 0-direction. Values of NC(LREG) and MC(LREG) 
are input into the program. Within the present DIMENSION statements, 
the maximum value of MC(LREG) is 19 and the maximum value of 

NC(LREG) is 40 with the additional constraint of 

8 

NC(LREG) = 150 . 

LREG=1 

Thus, the maximum number of mesh points used in the nacelle computa- 
tion is 2850. This does not include the points within the plume 
which will be discussed separately. Region 8 corresponds to the nose 
region. 

Increasing the number of mesh points in each region increases 
the resolution. To simplify the choice we have built into the pro- 
gram two sets of mesh distributions: one set for the complete nacelle 

calculation and the other set for the boattail option. The values of 
NC and MC of any. region may be input to change any of the preset 
values. The values of NC and MC for the nacelle calculation are 
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For the boattail 


NC(1) « 

18 

MC(1) « 

18 

NC (2) = 

18 

MC (2) = 

18 

NC (3) = 

25 

MC (3) = 

18 

NC (4) « 

16 

MC (4) = 

18 

NC (6) - 

18 

MC (6) = 

18 

NC(7) = 

32 

MC (7) = 

18 

NC (8) = 

19 

MC (8) = 

9 

option 




NC (1) = 

20 

MC (1) = 

15 

NC (2) = 

1 

MC (2) = 

1 

NC (3) = 

1 

MC (3) = 

1 

NC (4) = 

1 

MC (4) = 

1 

NC (6) = 

20 

MC (6) = 

15 

NC (7) = 

20 

MC (7) = 

15 

NC (8) = 

1 

MC (8) = 

1 


Note that all the values of NC and MC are not independent. 

) 

The constraints are 

MC (4) = MC (3) - MC (2) 

MC ( 7 ) - MC (6) = MC (1) 

NC(2) - NC(1) 


Once the number of mesh points are selected, the actual mesh point 
distribution in the physical plane is determined from the values 
of the stretching parameters. The stretching parameters x q , x-^, 
^2 > y Q , y 2 > a i> a 2» a 3 ’ an< ^ a 4 are determined from physical 

mesh point locations DD(J) , J = 1 -+ 6. The lengths DD are shown 
in Fig. 3.5. Shown here are the regional interfaces and the first 
and last grid lines (in the physical space). Values of DD are 
as shown in Fig. 3.5. Again, to simplify the running of 
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Fig. 3.5 Stretching Parameters 



Fig. 3.6 Plume Geometry 
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the program, two sets of values of DD are built into the program. 

* 

Complete Nacelle Calculation Boattail Option 


DD(1) 

« 0.005 

DD(4) -4.0 

DD(1) - 0 

DD(4) = 10 

DD(2) 

- 0.005 

DD(5) - 0.005 

DD(2) - 0.05 

DD(5) - 0.05 

DD(3) 

- 0.005 

DD(6) - 4.0 

DD(3) =0.05 

DD(6) = 10 


Note that the numbers 1-6 are only an index and do not correspond 
to any particular region. 

3. 1.1. 4 Plume Geometry 


Initially, the shape of the plume in Cartesian coordinates 
is input into the program. These coordinates, and y , are 

scaled to the initial radius of the jet as in Fig. 3.6. A piece- 
wise parabolic fit is made to the coordinates extending to down- 
stream infinity. After the plume computation has been made, the 
new boundary shape is curve fit automatically. 


3 -1.1.5 Flow Field Output 

There are two basic output options, the surface flow field 
output and the complete nacelle flow field. For the surface flow 
field, values of the Cartesian coordinate position, x,y, pres- 
sure p, Mach number M, flow deflection t = v/u, and pressure 
coefficient Cp are given for each mesh point along the surfaces 
of the outer cowl, inner cowl, and plume boundary (outside the 
plume) I' ' For the complete flow field, output values of the entire 
field are given for x, y, p, a (speed of sound) , u, v (velocity 
components), M, p (density), C and S (entropy). The output is 
divided into regions as in Fig. 2.2. Note that for regions 1-7, 
u,v are the Cartesian velocities and for region 8, u,v are the 
polar velocities. All the outputs are nondimens ional as discussed 
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in Section 2.2. The pressure is scaled by p^, the density by 
p , all velocities by V p^/ p^' , and all lengths by L, the 
length of the nacelle. 

In addition, the mass flow going through the inlet is numeri- 
cally integrated (trapezoidal rule) and output whenever the com- 
plete flow field output is called. 

3. 1.1. 6 Tape Input and Output 

The program has a restart capability. At the end of a given 
run* the contents of all common blocks are dumped onto Tape Unit 12 . 
At the start of the next run, the data are read from Tape Unit 11. 
The implementation of input parameters to use this option is dis- 
cussed in Section 3.1.2. 

3. 1.1. 7 Plume Calculation 

The plume is computed every KPLUME number of time steps . If 
KPLUME = - 1 there will be no plume and a straight pipe will be 
assumed to extend downstream of the body. The plume routines have 
been written as a separate independent program. It is made part 
of the entire nacelle program through subroutines BOUND and PLUBO. 
BOUND takes the values of the pressure along the plume outer sur- 
face (obtained in the nacelle calculation) and converts the re- 
sults into a form suitable for the PLUME routine. PLUBO curve 
fits the computed plume shape into a form compatible with the 
nacelle calculation. The plume is computed to a specific value 
of length which is input. This length should be shorter than the 
distance to the Mach disk. 

The plume routine output values of Z, r, p, u, w, 8 , M, S, 
and t, where quantities are nondimensional with respect to nozzle 
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exhaust plane conditions (Section 2.2.4). The major input condi- 
tions are M, y , PRATIO. Note that here PRATIO is the ratio of 
static pressures of the jet to the free stream. This quantity is 
related to the stagnation pressure ratio P f by 


PRATIO = 




{a . ^ »V' V ' 

|(1 + M 2 )^ 


1 . 

J 00 

1 . 

jet 


[Note that the stagnation temperature ratio is only required to 
rescale the variables to free stream conditions as in Eqs . (12)]. 


3. 1.1. 8 Over-all Logic Flow 

The details of computer program logic are presented in 
Sections 3.2.1 and 3.2.2. Here we will describe some of the im- 
portant features of the logic flow. 

The program begins by receiving input data on the basic 
parameters associated with the run, all options, and mesh point 
distributions. Following this, initial data for all variables are 
computed. The coordinate stretchings are found in subroutine 
STRECH. The values ot the body geometry are calculated in WALL 
and the initial value of the plume shape is given in PLUBO. The 
interpolations for the nose region are set up in SETNOS . 

Then the major time loop is entered. The size of the new 
time step is determined and the value of time incremented. Within 
the time loop, the entire calculation for the cusped nacelle is 
performed in POINT. A time sub-loop is set up for the nose region 
and is computed in NOSE. If the option for the plume is in effect 
and this is the appropriate time step for its computation, sub- 
routine PLUME is called. This routine acts as a main routine for 
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the plume computation with the main calculation taking place in 
SUPER and the plume output given in OUTP. Once a new plume shape 
is determined, it is curve fit in PLUBO and the stretchings for 
the nacelle calculation redetermined in STRECH. Now, except for 
some optional output for the nacelle in OUTPUN, the step of the 
major time loop is completed. 

After a specified number of time steps, KA , the computation, 
is completed. Finally, there are optional outputs in OUTPUN and 
optional saving of the data on tape in TAPER. 

3.1.2 Program Usage and Operation 

In order to explain the usage of the program, we present a 
summary and description of all input data. The order of the dis- 
cussion follows the order of the data cards. Examples of all data 
input are shown in Section 3.2.4. Some of the input is of the 
NAMELIST type. This input procedure has the advantages of being 
easily identifiable and does not require all the data mentioned 
in the namelist declaration to be input. Thus, much of the input 
is optional and these paraemeters will be denoted by an asterisk 
(*) . The values of the optional input parameters are denoted de- 
fault options and will be discussed below. 

3.1.2. 1 Program Input 

NAMELIST/RUN/NRUN, NDATE, EM, GAMMA, RMFLO, LA, LSYM, KPLUME , STAB 


NRUN 


Run number, Integer 1-99,999 

NDATE 


Rundate, dimensioned NDATE(3), typical input 



NDATE =2,19,73 

EM 

M 

00 

Free stream Mach number 
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GAMMA* 


Free stream ratio of specific heats, default 
value 1.4 

Mass flow ratio at engine-inlet interface 

= 0 two dimensional, = 1 axisymmetric 

= 0 complete nacelle, = 1 boattail option 

The number of time steps between plume computa- 
tions, typically 100. If equal to -1 no 
plume computation will be made, default value -1 

Stability parameter (for t ime step determination) . 
Typically 0.75-1.0, default value 1.0 

3. 1.2. 2 Input/Output Parameters 

NAMELIST/ INOUT/MRE AD, MRITE , KA, JA, JB, MB, LOUTl, LOUT 2 

Tape read parameter, integer value represents 
number of read cycles, =0 no tape input, 
default value 0 

Tape write parameter, = 0 do not write on tape, 

= 1 write on tape, = -1 rewind tape and write 
on tape, default value 0 

Total number of time steps 

Number of time steps between surface flow, 
field outputs, = -1 for no output 

JB* Number of time steps between complete flow field 

outputs, = -1 for no output, default value -1 

MB* Number of mesh points from nacelle surface in 

complete flow field output, = -1 all mesh points 
default value = 1 (surface values only) 

LOUTl* =1 for surface output at end of computation, 

= 0 no output, default value 0 

L0UT2* = 1 for complete flow field output at end of 

computation, = 0 no output, default value 0 


MREAD* 

MRITE* 

KA 

JA 


y 

RMFLO m 

r 

LA 

LSYM 

KPLUME* 

STAB* 
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3. 1.2. 3 Mesh Point Parameters 


NAMELIST/MESH/NC, MC, DD 


NC* 

MC* 

DD* 


Number of mesh points in x-directions , dimen- 
sioned NC(8), full discussion in Section 3.1.1 


Number of mesh points in y-directions , dimen- 
sioned MC(8), full discussion in Section 3.1.1 

Mesh point distribution parameters, dimensioned 
DD(6), full discussion in Section 3.1.1 


3. 1.2. 4 Geometry Input 

1. General Description A Format (20A4) 

A Alphanumeric description, 80 characters long 

2. General Parameters ELL, HBAR, RNBAR, JNOSB, JOUTB, JINB 

Format (3F10.4, 3110) 

ELL L Length of the nacelle 

HBAR H Height of nominal center of cowl_lip above 

centerline, same dimensions as L 

RNBAR r Radius of cowl lip (nominal), same dimensions 

1 as Z 

JNOSB Number of divisional points to be input for 

cowl lip, = 0 for boattail option 

JOUTB Number of divisional points to be input for 

external cowl 

JINB Number of divisional points to be input for 

internal cowl, = 0 for boattail option 
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3. 


Cowl Lip Data XBAR, YBAR, YPBAR Format (3F10.4, 3110) 


- 'v 


XBAR 

X 



YBAR 

y 

. Cartesian position and slope 
points along cowl lip. Same 

of divisional 
dimensions as 

YPBAR 

dX 
dx j 



There 

will be 

JNOSB input cards of this type. 

External Cowl Data XBAR, YBAR, YPBAR, LFUNO(J) 
Format (3F10.4, 3110) 

XBAR 

X 



YBAR 

y 

, Cartesian position and slope 
r points along external cowl. 

of divisional 

YPBAR 

dZ 

dx j 




3 for cubic fit between J and J + 1; 

1 for straight line between J and J + 1 


LFUNO(J) 

There will be JOUTB input cards of this type. 


Internal Cowl Data XBAR, YBAR, YPBAR, LFUNI(J) 
Format (3F10.4, 3110) 


XBAR x 


YBAR 

YPBAR 


Cartesian position and slope of divisional 
points along internal cowl. 


& 

dx 


LFUNI(J) 


= 3 for cubic fit between J and J + 1; 

=1 for straight line between J and J + 1 


There will be JINB input cards of this type. 
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3. 1.2. 5 Plume Geometry Input (Initial) 


1. General Description AW Format (20A4) 

AW Alphanumeric description, 80 characters long 

2. Number of Intervals KK Format (1015) 

KK Number of intervals for plume geometry. Corre- 

sponds to the number of input cards to follow. 




If KK = 1 the plume will be a straight pipe 

Plume 

Shape 

XPP, YPP Format (2F10.4) 

XPP 

X 

p 

Cartesian coordinates of initial plume shape. 
Scaling shown of Fig. 3.6. There will be KK 

YPP 

y»'l 

of these cards . 


3. 1,2. 6 Plume Calculation Input 


NAMELIST/ JET/NA, MA, KA, KOUT, PRATIO, DIST, GAMMA, STAB, ACH 


NA* 

KA 

KOUT 

PRATIO 

DIST 

GAMMA* 

STAB* 

ACH 


Y 


M. „ 
jet 


Number of mesh intervals in the r-direction 
NA(1) low pressure side of the shock, NA(2) 
high pressure side of shock, maximum value = 30 

Maximum number of AZ steps in calculation, 
not crucial, calculation should stop at Z = DIST 

Number of steps between plume output 

Static pressure ratio p, _/p . Discussed in 
Section 3.1.1 

Distance to end of jet computation. Should be 
before Mach disk forms 

Ratio of specific heats for jet 

C-F-L parameter, usually = 1 

Mach number at exhaust plane of nozzle 


A schematic of the card data input is shown in Fig. 3.7- 
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NAMELIST RUN 













3.1.3 Accuracy and Limitations 


The basic numerical methods for this computation are essen- 
tially second-order accurate. However, the boundary conditions 
and the interpolations at the nose region tend to reduce this level 
of accuracy. Another factor affecting the accuracy is the spatial 
mesh point resolution. The resolution is determined by the number 
of mesh points and the stretching parameters. Improving the reso- 
lution by increasing the number of mesh points increases the com- 
putational time and core storage. The values of the stretching 
parameters suggested as default options in the computer program 
should be sufficient for most applications. In general, the only 
effective means of evaluating the accuracy of a computational ap- 
proach of this complexity is by comparing with other analytical 
methods and with experimental data. Caution must be exercised in 
comparing an inviscid computation with data because of viscous 
effects. Regarding the nacelle computation, an indication of the 
accuracy of the approach is given by the data comparisons in 
Section 2.3. ' 

There are several limitations to this computational program. 
Namely, accurate results can only be assured for subcritical free 
stream Mach numbers. Furthermore, even at high ' subcritical Mach 
numbers, when the mass flow ratio is low, problems may develop. 

The reasons for these limitations are discussed in Section 2. Since 
there are no provisions for a bow shock, the program will not work 
for supersonic free streams. 

The basic program is written to handle a cusped centerbody. 
However, the geometry routine would have to be modified to run 
this case. In addition, the program cannot handle the short cowl 
nacelle. An approximation of the short cowl can be made by treat- 
ing the plume with a specified internal solid boundary which 
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represents the jet plume. The remainder of the plume would corre- 
spond to the fan jet. This computation requires a change in the 
plume geometry input procedure and the plume computation routine. 

3.2 PROGRAM ORIENTED DOCUMENTATION 

In this section, we attempt to present enough detailed in- 
formation about the computer program to enable the user to under- 
stand and possibly to change the source language code. Firstly, 
the structure of the program is schematically described by flow 
charts of the main program and all major subroutines. A sub- 
routine tree diagram is also presented. Then each subroutine is 
discussed followed by a description of input/output files. Lastly, 
the input and a partial output for two test cases are presented 
along with a FORTRAN source program listing. 

3.2.1 Program Flow Charts 

The flow charts for Program 15C are shown in Figs. 3.8-3.18. 

A subroutine tree diagram is presented in Fig; 3.19. 

3.2.2 Subroutine Description 

POINT Computes one time step for all interior points in 

the cusped-nacelle formulation of the problem. All 
solid boundary conditions, free stream conditions, 
and interface matchings are handled here. Included 
is a separate computation for the cusped-nacelle 
point. 

STRECH First calculates all coordinate stretching parameters. 

Develops coordinates and derivatives of mesh points 
for cusped nacelle formulation in the physical plane. 
Obtains information from WALL and PLUBO. After major 
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Fig. 3.8 Main Program 
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EQUATIONS 
PT, UT, VT. ST 


INTERFACE 

' 

INTERFACE 

MATCHING 

(FINAL) 


< 

MATCHIIMlj 

(INTER.) 


f 


t 

FREE STREAM 


FREE STREAM 


VALUES 


VALUES 


(INTER.) 


(FINAL) 




CUSP POINT 
COMPUTATION 


I 


END OF LOOPS 
ON LREG, N. M 







BOUNDARY 

CONDITIONS 

(INTER.) 


BOUNDARY 

CONDITIONS 

(FINAL) 

i 


' i 

i 

INTERMEDIATE . 

VALUES 

PI, Ul, VLSI 


FINAL 
VALUES 
P, U, V, S 

: 3 


1 

k 




Fig. 3.9 Subroutine Point 
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Fib. 3.10 Subroutine Strech (MREG) 











INPUT 
DESCRIPTION 
AND GEN. 
PARAM- 
ETERS 



COMPUTE 
COWL LIP 
PARAMETERS 



INPUT 
INTERNAL 
COWL 
DATA 


COMPUTE 
INTERNAL COWL 
PARAMETERS 



INPUT 
INTERNAL 
COWL 
DATA 


COMPUTE 
INTERNAL COWL 
PARAMETERS 



CENTERBODY 

GEOMETRY 


I 


RETURN 
Y, YP 


COWL LIP 
GEOMETRY 


I 


RETURN 
Y, YP 


3 


INTERNAL 

NACELLE 

GEOMETRY 


I 


RETURN 
Y, YP 


EXTERNAL 

NACELLE 

GEOMETRY 


RETURN 
Y, YP 


Fig. 3.11 Subroutine Wall (I REG, X, Y, YP, LP) 
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Fig. 3.12 Subroutine Plubo (x, y, yp) 
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ENTER 

JRITE 


SET UP INTER- 
POLATION 
PARAMETERS 
r = c 



NO 


OUTPUT 
^PARAMETERS/ 


SET UP 

INTERPOLATION 
AT INTERIOR 
POINTS 



NO 


NO 


OUTPUT 
^PARAMETERS/ 


OUTPUT 
k PARAMETERS/ 


Fig. 3.13 Subroutine Setnos (JRITE) 
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C ENTER IREG. 
NTIMES.ITYPE 



YES 


YES 


COMPLETE 
FLOW FIELD 
OUTPUT 

; 

f 

INLET 

FLOW 

INTEGF 

AASS 

ATION 




NACELLE 

SURFACE 

OUTPUT 


■ 1 

r 

GEOM 

OUTPL 

ETRY 

JT 



^ RETURN 




Fig. 3.15 Subroutine Outpun (IREG, NTIMES. ITYPE) 
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ENTER 

IREAD, IRITE 



Fig. 3.16 Subroutine Taper (I Read, I Rite) 






Fig. 3.17 Subroutine Plume 
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ENTER 



CALL 

BOUND 


Fig. 3.18 Subroutine Super 










STRECH. 

SETNOS 

oUtpun 

TAPER 


PLUBO 


PLUME 


BOUND 
OUTP 
SUPER - 


POINT 


NOSE 


Fig. 3.19 Subroutine Tree Diagram 


BOUND 
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computation begins, this routine modifies the co- 
ordinates after each new plume computation. 


WALL 

PLUBO 

SETNOS 

NOSE 

OUT PUN 

TAPER 

PLUME 

SUPER 


BOUND 


OUTP 


Performs curve fit to nacelle geometry. 

Performs curve fit to plume shape, either from in- 
put data or as the result of the plume computation. 

Sets up all interpolations used in the nose region 
computation. 

Computes all points in the nose region. Performs 
all small time steps for one complete time step. 
Handles all interpolations at the boundaries . 

Handles all program output from the nacelle calcula- 
tion (excluding the plume) . 

All tape input/output for restart capability. 

Complete plume computation MAIN routine . 

Computes all points for the plume calculations. 
Handles shock points, axis points and plume boundary 
points. Predicts and fits the imbedded shock. 

Takes static pressure solution along plume boundary 
from nacelle calculation and converts to plume non- 
dimensionalization . Interpolates pressures for 
plume boundary computation in SUPER. 

Handles all output from plume computation. 


3.2.3 Input/Output Files 

The entire program can be run using the standard input/output 
files (Tape 5 for input and Tape 6 for outputs) . During the com- 
puter program run, no intermediate tapes or disks are used. However, 
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the restart capability of the program does utilize tapes or perma- 
nent disk files . 

All output to the tape or disk is on TAPE Unit 12 and all 
input on TAPE Unit 11. Schematically, a typical sequence of runs 
may be as follows 


NRUN 

MRE AD 

MRITE 

INPUT 

OUTPUT 

1 

0 

1 

TAPE 5 

^TAPE 12 

2 

1 

1 

TAPE 11^ 

^/TAPE 12 

3 

1 

1 

TAPE 11^ 

^/TAPE 12 

4 

1 

0 

TAPE 11^ 

- 


Note: MRE AD and MRITE are parameters for subroutine TAPER and are 

described in Section 3.1.2. Also standard output on TAPE 
Unit 6 will be produced for all the above runs as discussed 
in Section 3.1. 


3.2.4 Test Cases 


Sample Input - Nacelle Calculation 

SRUN NRUN=1 12*NDATE=1 » 8 » 7 3 » EM= • 9 ,RMFLO= .8852 * L A= 1 .LSYM=0 SEND 
SINOUT KA=400 , J A= 1 0 0 . JB= l 00 , MB = -1 .LOUT 1= t «LOUT2= 1 ,MREAD=0. MR I TE = 0 


smESH 

S 



I NLET 

NO. 8 

NACA 1-85-100 (C. 

R.=l .093) 

18. 

7.682 

0.200 

4 4 

"■ .200 

7.825 

~ “ o;3Y?rs7 


.072 

7.770 

.6111 


: 0 • o 

, 7,682 

1.0 • , 

, 

.~25o 

7.455 

-.4270 


.200 

7.825 

.37037 

3 

.7 20 

7.966 

0.2222 

3 

■ 2.? 

8.280 

0 • 1 22~2 

3 -- - 

■ 18. 

9. 

0 . 

3 

.200 

7.455 

— i4270 

3 

” 0TT5T~ 

7TT48 

0.0174" ‘ 

1 "• 

4.5 

7.4 13 

0.0174 

3 

8.10 

7.527 

0.0 

3 

— NO 

PLUME' 
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Sample Input - Boattail/Plume Calculation 


$RUN NRUN=108,NDATE=1 , 22 , 7 3, E M= . 7 , L A= 1 ,LSYM=1 ,KPLUME=50$ 

, $INOUT KA=500.JA=-1,JB=50,MB=4 $ 

6 )^2'0T^CT'7 )^20^wcm =1 ~ */ 

BO AT t AIL. STRAIGHT PIPE 

1 • • 5 ' 0. 0 1 0 

xr cr. ~~ i . — 

PLUME DATA 

T :'2 ' ' j ■ ' ■ 

2Y T.T5 ~ 

4. 1.24 

$JET KOUT= 1 0. PR AT I 0=3. ,TTOT=l,KMAP=l 0, KA = 1 90 ,NA< 1 ) = 1 5 , NA (2) = 10*DIST = 4 < 


$ 


sF.Nr. 
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Sample Output - Nacelle Calculation 
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3.2.5 Source Language Listing 


_P-RCGRAM..NACtL4INPUI.nUTPUTtTAPE5=IiuPUI^.TAPE6aOUTPUT-»TAPEl 1 ,TAPEl2> 
TRflNSONIC NACELLE CALCULATION 
THlKt,8LUNT Cowl L^P 

_COM. m CN/HLN.1/’NC (8> «KC'(8V ♦NCI »NC2»NC-**-NC4*NC5*NC6»NC-7*NC8»MCi *MC2 


1 ♦MC3fMCA < MCPfMC6fHC7,McA*NfiEG (8) ,NiuC (2) ,Mmc (80 • 2 ) *NMAX*MrAx 
2,GAMM*GA.GB t GC*G0,Ge-*GF,X <4o,8) , Y , 1 9 ♦ 8 ) ♦ XXP ( 1 30 ) 

-3-»XVP < 130*1? > *HH*xE.Yr* YA*XC*R o»RO*oMFLO»TT*CC*EM»PH~ — 


4, SI* THE (20) *COSTHE (2* ) . R ( 20 « 1 9) * IS vM *L A *D X < 8 ) *0Y (8) 
COmMcN/GL^S/P (ISO* 1?', ♦U(150*19) * V ( 1 50 ♦ 1 9 > »S<150»1?) *P1 (150*19) 
_l.iUUiS'4» 1.9.1 .*V.I tlSO.lo) *SL< 15 q *19 )*/yS(d »2)*NF(8 *2 ) 


2*MS<8*2)*MF<8,2),TIMr*OT,K*J,XOTT(?) »Q INF . 0 INFN *KDT VS * OTS 
C08MCN/8LN3/YET ( 130 * T 9 ) , XCS ( 1 30 ) . Xo * X 1 , X2 , X3 ♦ YO ♦ Y2 ♦ ALP ( 4 ) ♦ DO < 6 ) 

4 * BET (4 1 *lSLP*HUM3o)'.HL( 13Q) *HCn3«) »HiiPR(13o» fPLPP(130^^HCPR(130)NAC00l20 


NaCOOOIO 
NAC0O020 
NaCOoo30 
NAC01040 
NAC00050 
-NACOoOftO 
NaCO 0070 

nacoooso 

-NAC00090 
NACOO 1 00 
NACOol 10 


COMMCN/BLK<t/EPl (20 ) *fP 2<20) *EP3(10*> «EP4 < 1 00) *EP5 «38> *EPfc <38) 
l*NNl(2j)*N\<i(20)*Hl(^0)*P2(20)*M3(T8),M4(38)*Ll(lOO)*L3(100) 
-Zt-llAlOalilLLlOO)-. Jfi.(a0,J9.MB.(20UBDR,(2A4_ — - 1 


C0PfCN/BLK6/ XPL ( 2 0 0 > *PPL (2 00) tKPP 
COMMcN/dL^P/ PRAt1C*oOIST*PRAD*KPLiiME* JJ 
-CGM^CNt/BL^R^OfVtWwtSTABtNNMAXtHTFfrNMTh 


DIMENSION D ] T (7 ) *NOArE (3) 

NApt L IST/RLN/ npUN*NOaTEiEM*GAMMA,RmELO.LA.lSYM, KPLUME.STAB 
-NAwEl IST /_I\yu.T/MREAO M r I TE *KA*-JA-* Jo *M8*-L0wT4 *LOL ; T2 


NACO 0130 
N ACO A 1 40 
-NAC00150- 
NAC00160 
NACOO 170 

-naCooibo - 

NaCOO 190 
NAC00200 
-NAC00210- 
NaCO 0220 
NAC00230 


NAKELlST/MESH/^CtMCfnO 
CALL ERRSET (208*256. .1,1 ,0*0) 

JM104— FORMAT (1H1 *29X *1 1 HPRaGrAM_JLSC//2oX,29HTRAnSON I C— NACELLE COmPut ATIONAC00240 
IN //20X .3MRUN* 15, 13 X, 12, 1 M/* 12* In/, I 2//20x-r5HMACHc.F5. 3 *eX,5HPFLO=NAC00250 
2.F6.4///4X, 1 OH INPUT nATA/4X t2HNC,5», 014/4 x f 2HPC*5X, 8 14/4 X»2H0D. NAC00260 

3-*-5X-»6F8^3/4 i X,6HMREA0 ! l.-I2»3X.»6HMRI.Tps-*4?*-3x*-3HKA*^J5*3X*3pJAs t I3 — NaC00270 

4,3x*3HJB=, I3/4X,7HKP; UME=»14/1 NaC0O?80 

1002 F0RMAT(4X*2lHTW0-DIMrNsi0NAL FLOW /) NAC00290 

FORMaT (4X« 1 /HAXIsYKMcTRIC FLOW-/-) : NaCOOSOO 

F0RM.aT( 4X» IbHPoATTAIj OPTION /) NaC0o310 

FORMAT (2X*Ib*2X,7 (I2',E10.3) ) NAC00320 

JNPyT.JDATA __ , NACO 0 330 

G AP.M A s l » 4 NAC0ri340 

R M F L o = 1 • NAC0O350 

STABa 1 . = _ — : — — NAC00360 

KPlUwE=-1 NAC0O370 

REA^(5*RUN) NaCo 0380 

mbemdso — r -NaCo 6 390 


-1003 

1 004 

1 005 


MRITe= 0 
MBs 1 
JB="1- _ 


LOUT 1 =0 
LOU l ?=Q 

_REAP(5,IN0UTL 


IF (MPEAO.NE.O) GO TO 80 
IF (LRYM.EU.I) GO TC TO 

_MC ( 1 ) a 1 8 : 

MC(2)=)8 
NC ( 1 ) = 1 8 

_NC (3) =25 

N c ( 4 ) a 1 b 

NC(6)al8 

~NC ( 7-i-«32- 


NAC00400 

NaCo 04 1 o 

-NAC00420 

NaCo 0430 

NAC00440 

-NAC00450 

NAC0O460 

NAC00470 

-NAC00480 

NAC00490 

NAC0O500 

NACOoSlO 

NAC0O520 

NAC00530 

-NAC00540 


100 



NC< a , » 1 9 

MC<8)=9 
00 ( 1 ) =.005 

•OD(2)».uOS ; 

00<3) = ,c'05 
00(4)s4. 

■00-4 5) c , OOS - — — — 

0D(6)=4, 

60 To 20 

i >=lv 

mc ( 2) s] 

MC<8)=1 

-NC 4 a > = 1 0 : 

■ NC ( 3 j a 1 . 

NC < 3 ) =1 

— — 

NC (6) =1 j 
NC ( 7 ) a ) 0 

nc-< a ).= i : ; 

oo ( i j =o . 

00 ( 2 ) s , '>5 

-00 (3.) a,;: 5 

0 0 < * ) = l 0 . 

00(5)3,05 

. 00.(6) si J._; _ __ 

2c REAU (5. MESH) 

JJcLA 

CONSTANTS _ 

MC(3)=MC<2> 

MC(*)=MC(2) 

.MC (5) aM.C (.2.) 

MC (6 ) = MC { 1 ) 

MC ( 7 ) aMC ( 1 ) 

NC (2) aNC-lU 

NC (5) a) 

WRITE (6*1 O ') 1 ) NRUN,NnATE»EMtRMFLO*K'C*MC»DD»MREAC*MRlTE»KAf JA*J 8 
L 4KP.L1.ME 

IF(La.EQ.O) WR'ITE (6*7002) 

IF(La.EO.I) WRlTf (6»7003) 

IF. A L S Y. M • E.H 4 ..Ll W SJ I E I *j. 1 00-4.) 

NCl=NC ( 1 ) 

NC2=NC (2) 

NC3=NC(3) 

NC4=nC (4 ) 

NC5=NC (5) 

'NC6?NC ( 6.1. 

NC7=NC(7) 

. NC8=NC(B) 

— —MC I = MC ( 1-) 

MC?=MC (2) 

MC3=MC(3) 

MC4=MC <4 ) 

MC5=MC(5) 

MC6=MC(fa) 

MC7 = MC 17) 


TOT 


NAC0O55O 

NAC00560 

NACO-1570 

NACO 1580 

NAC01590 
N a£ 0 160 0 

NAC01610 

NAC00620 

NAC01630 

NAC01640 

NACOlfcSO 
NAC0O680 
— — NAC01670 
NAC00680 
NaCO 1690 ’ 
- — NAC01700 
NAC00710 
NAC00720 
— NAC00730 
NAC007A0 
NAC00750 
— NAC01760 
NAC01770 
NACOnTSO 
— -NAC01790 
NaC‘ 01800 
N ACO i 8 1 0 
- — NAC0i820 
NACO.1830 
N A CO i p 4 0 
— NAC00850 
NAC00860 
NAC01870 

NAC0O880 

NACO 1890 
NAC01900 
NAC0191O 
NAC0O920 
NAC 01930 

NAC00940 

NAC01950 

NAC01960 

NAC0O970 

NAC01980 

NAC01990 

NACOlOOO 

nacoioio 
NACO l 020 
— NaCOI 030 
NACOIOAO 
NAC01050 

NAC 01060 

NACO 1070 
NAC01080 
NAC01090 



HC{8 ) 

(1 ) =0 

0 LREG = <J»0 

<LR£GJ sNREGlLREG-D ^NtaREG-l )- 
.11=2 
*21=2 

* U =2 ^ 

• 21=2 
* 1 ) =2 

• 2 ) =1 

.11=2 
* 2 > = 1 

♦ 1 ) -2 

.21=1 

. 1 ) =2 

.2 1 = 1 

1 1 ) =2 
« 2 ) = 1 

• I 

1 2 ) =2 

• l ) =nci 

. 2 i=nc i-i : 

1 1 ) = NC2 

.2) =NC.3-1 

«.u=Nca 

» 2 ) = N C 3 - 1 
« 1)=NC4 

.21-SNC4 

, \ ) SNC5-1 

• 2 ) snCS — 1 

» 1 )=nC6-1 

,2j=NCe-i 
1 1 ) =NC7 

» 2l=NC7—l : : 

1 1 ) =NCP 
.2) =NCP-1 

i U-=l 

.21=1 

.11=1 

.21=1 

.11=1 

.21=1 

. I ) = 1 

.21=1 

.11=1 

.21=1 

.11=1 

.21=1 

. 1-) =1 

. 21=1 

.11=1 

.21=1 : : 

, 1 1 =mC 1 - 1 
1 2 1 =mC 1-1 

. d=mc2 — h>i 



*c 

c 


c 


M F-< 2 , 2 > = M C ?-l NACo 1 650 

MF(3,1)=mC3 NaCO 1 660 

MF (3,2) =MC3 NACO 1670 


+4F (4, l >=MC4 — : NACO 1680 

MF(4,2)=MC4 NaCO 1 690 


MF <5, 1 ) =mC' 

-MF <5,2> =MC5 — 
MF fb, 1 ) =MCc-l 
MF (b,2) =MC6“1 

-WF-(7,1 )=MC7-1- 


MF ( 7 ,2 J =MC7-1 
MF (8, 1 ) sMC» 

-MF (8,2) iMCc-1- 

««<•«««««««<«» 
the following three 


NAC01700 

NaCO 1710 

NACO 1720 
NAC01730 

NACO 1740 

NACO 1 750 
NACO 1 760 

N ACO 17 70 

NAC01780 

rAROS MUST BE CHANGED ACCORCING TO DIMENSIONS NAC01790 


-NMA A*40~ — 

NNM«X=1 50 
MMA*=19 

a« «*««««•«« * a « 

GAsGaMMA/ <GAMMA-1 . ) 
G8=1./(GAMva-i.) 

_GC=.(G AMM.A +.1 . 1 /< GAMMA. l^T 
60= (GAMMA- 10/2, 

GE = <GAMMA+ 1 . ) /2. 

-GF = Si;RT (GAMMA) 

QO=Gf«Em 


NACoiROO 

NaCOIRIO 
NaCO 1820 
— —NACO 1 R30 
NAC01840 
NaCO 1 850 

NaCO I860 

NACO l R70 
NAC0IR80 

NAC0I890 

NACOIROO 


PTI a 4.*ATAN (1.) 

jua 

j=o 

time =6. 

DJ.a l 

00 4p LREG=ltH 
IF(Nc(LREG) «LE . 1 ) GO To 35 

IF_(Mc (LHEG) .LE.a ) GC -T0_3S 

OX (LREG) =1 ./FLOAT <NC <LREG) -1 > 

oy (Lregi =1 ./float (mc',lreg) - i j 

GO-J.C-4 0 

35 OX (LfiEG) =0. 

OY (LpEG) =0. 

AO CONI I NUE 

DO 50 LREG=T *8 
00 4j m=1»mMax 

4_2_ Y. ( M.j L P E G ). = .( M l-l ) *0 Y X L o E G J 


NACO 1910 

NACO 1 920 

NaCO 1 930 
NaCO 1 940 

NAC01950 

NACO 1960 
NaCO 1970 

NaCO 1 9R0 

NACO \990 
NAC02000 

NAC02010 

NAC02020 

NAC02030 

NAC02040 

NAC02050 
NAC02060 
— NAC02070 


00 4 A N=1 »\MAX 
44 X (N»lREG) = (N-i ) * 0 x (L oEG) 

5X,_C0nT INUE_ 

LSlt = 0 

NTh=NC4/2 

NNTMsNREG (A ).»NTH 

KPP=kC6 

CALh WALL H *0. * YO«OUwtl ) 

YA. =10 

CALL WALL ( 6 ’ 1 *. . ♦ yD . DmM * 1 ) 

YE = YO 

XO r AE-— 


TUT 


NaCo?o80 

NaC02090 

NAC02100 

NAC021 10 
NaC02 1 20 

N AC02 1 30 

N A C 0 2 140 
N AC02 1 50 

NAC02160 - 

NAC021 70 
NAC02180 
NAC02190 



CAUL WALL L7 »XD#YDjCUmJ 

PRaUsYD 

PDlSTaxE 

IALL.T I ALI2A J.lQM : 

00 60 Ns 1 t NNm A X 
00 6q Msl,^MflX 

P.lNtM ) =0. 

U(N«MJ=(,. 

V (N»M) *C . 

S-( N»R1 = U . : 

PI (NfM)so. 

UI (N.M) =0. 

Vi (N»M) =0 • 

C SI(N»M)=0. 

NN4=NREG (4) *NC4 

00_I(i_M= U-RtA 

P INN* ,M) s-2 .3 
5 S(nN4»m)=0. 

— WWs *5 

OlN^aO. 

Q I N E N = 0 • 

— *&^UU=o. 

X0Tt(2)s 0 . 

— CAtL—STWECHiSi 

IP (LcYM.NE.i) CAUL ScTNOS (L0UT1 ) 
IF(LoUTl.NE.O) CALL nUTPUN < 1 » 1 • 3 ) 


GO To 100 

call taper ( m reao,mriH) 

_.M A I tv -LOOP 

K = K ♦ 1 

time STEPSlZE DETERMINATION 

_C-F--L-RULE — 

0T=1 . 

00 120 LREG=1,7 

OTT ( L R E G > s 1 . 

IP (LREG.EO.b) GO TC T2o 
NCC=NC (LREC-) 

_KCC=MC (LREG1 

MREG=l 

IP (LrEG.E0.I.0R«L r EG’GE. 6) MREG = 2 
— IE-lLS YM , EO «-L« AND ,MWE<» «EQ»+-) — GO— TO— +20 

00 A l0 N=2»Ncc 

NNsNrPG (LREG) *n 

_IE_<MrEG .EQ.iiJ — 60_ JC--T-02— : 

M=MC2 
L = m 

Ga. To ic4 : 

Ms I 

L=M*l 

_A = gE«SOR-T (EXP (P (nN*mV/GA*S (NN.M) /-GaMMA>- 
OPA = SORT <u (NN, M ) 0*2 (NNfM) **2 ) ♦A 
OYRSABS ( YYF <NN,L1-YYr>(NN,L-l > ) 
_0XR=SQRI.UXAP_(NN)»XXt^U)N-lU-»«2M-YyP4NN 
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W I «J 


— OSsAmINI (L)xR,dyR> — 

0Tl=STAH*O=;/QPfl 
OTT(LREG)=AMlNi-(OTT (i'REB) »DT1 ) 

IF ( HP E G . E Q , 1 . A NO . H . En . MC 2+- 60-T0— IS* 

llo continue. 

120 DTs^mINI <DT *OTT llREG', ) 

I-F~(*V«LE • • 1 1 -OTbOT/2 % 

IF ( (K/1 .(■ ) *1 O.EO.K) WolTE (6>* 1005)K» (LREGtOTT <LREG) *LREG=1 ,7) 
IF{0T.GT.l.t-5) GO Tn 130' 

CALL 0i.iTPuVt9- r -l ,2) 

call exit 
l 3 q dxx=ox ( a) 

1 F-< <; s x M ^ E Q 60—TC — ?5o : 

DYY s DY<8) 

M=1 

OT S 1 » U : — 

00 140 N= 1 » NCR 
NNsNREG (H) * N 

A =gF « SO R T. ( E APT p < nN , m', /g A ♦S INN ♦ M MG * MHAT-l 

OP A = SORT <U (Nn,v ) «* 2*w (NN«MJ *«2) +A 
ROTHsR (N»M) *PI i«oXX 

DRRT-J CC-a ( a ) 1 *D Y.y- (Y_f M UtMLP I_I_Anx X 

DTl=STAb*AM!Ml < DRP , RnTH MQP A 
1*0 OTSUAHInI (U.TS1 .OTl) 

KO I Vc = 1 ♦ 0 T / U T-S l i 

OTS=OT/KOl VS 

IF < (K/.I-tj)*lO.EQ.K)h.*fTE(6«l0O5) k » kDI VS t DTS 

A Sul_T I M t= T I ME > C i- : 

C ACCELERATION FUNCTION 
LAC^sO 

I F_U I- M E . G E . 1.. / w W ).G C — t 0— t 6 p ; : 

0UH^wW*PH«Time 

XoTT (2)=-0. *WW«PlI«SrN(nUM)/2, 

, Q I N E A =.♦.<! 0 * . >*. LI . -CCS<DUM-U 

GO To 170 
16o LACisl 

X 0-T X < 2 ) . = 0. 

QInENsOL 
1 70 CONTINUE 

£ 

IF (LSUP.EO.*) GO TC T80 
LSUP=0 

MSUP= C - ; - 

1 7? NNT=NREG (4) *NC4 

GO Tc 176 

1 7.4_N N J. = N NT h ♦ 3 _ ; 

MSU p =l • 

176 AS=Gf*S(JRT (EXP (P (NNT'.MC4) /GA*S (NNT.MC4)/GAMMA ) ) 

OU« = SQRT <U<NNT.,MC4) *i2*V<NNT. t MC4y.«*2M-AS 

IF (DlH.GT * 1 *05) LSIP-LSUP*1 

if (esup.eo.O) go tc Tso 

IF ( i-S'JP « £0 . i • AND ,MSUn « EQ »-C- > — GO— TO— y-74 

IFfLsUP.EO.l. ANO.MSUo.EO.il GO TO lflO 
NC4 =N Th ♦ 3 

hC (4 ) =NC4 ; : 
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- NAC02750 
NAC02760 
NAC02770 

— NAC027R0 
NAC02790 
NAC02800 
NAC02810 
NAC02820 
NAC02B30 
— NAC02840 
NAC02850 
NAC02860 
— NAC02870 
NAC02880 
NAC02890 
— NAC02900 
NAC02910 
NAC02920 
-NAC02930 
NAC02940 
NAC02950 
— NAC02960 ■ 
NAC02970 
NAC02980 
-NAC02990 
NACO 3000 
NAC03010 
-NAC03020 
NaCO 3030 
NACO 3040 
-NAC03050 
NAC03060 
NAC03070 
- NaCO 3080 
NAC03090 
NAC03100 
-NAC03110 
N A C 0 3 1 2.0 
NACO 3 1 30 
— N A C 0 3 1 4 0 
NACO 3 1 50 
NAC03160 
-NAC03170 
NAC03180 
NAC03190 
-NAC03200-- 
NAC0321O 
NaCO 3 220 

- NAC03230 - 
NAC03240 
NACO 3250 

-NAC03260 
NAC03270 
NACO 3280 
-NAC03290 



♦ 


NF.X* «.L) =NC4 

NF <*.2> =NC4 
IbO COnT iNUt 


IF<KPLUME.‘-.U._ 1 ) GC }0 200 

IF( (K/KPLU''t)*KOLUvE.NE.K) GO TO 2n0 

DO~l90 N=1»NC6 

NNsNfiFG <6> *N 
PPL<N)sE.XP(P(NN.1) ) 

1-9T— XPL-lN ) =XX P-(-NNi 

CALL PLUMt 
CALL STHECH<6) 

20-0-C0N-T.INUE 

c 

CALL point 

1 F_lL S YM . N _c A LJ NaS E- 

c 

XOTT (1)=XGTT(2) • 

Q- Ifjt-aQlNFH— : 

C 

IF(JA.EQ.-l) GO TO 2?0 

IP ( lK/JA._).« JA-.EO.« >~CaLL— 0UTPUN<9.MB*-H 

2 1 0 IF(JH,EO.-n GO TO 2?0 

IF ( (K/JH) *o«.Fo.K) C/iLL OUT PUN (9 .Mn » 2 ) 

— 22n — IF (K.L T .Krt ) <jq tG— lOO — 

CAi.L T«PER(‘J.mRttE) 

IF (L0UT2.EG . 1 ) CALL aUTPUN ( 9 , MB « 2) 

STOP- — 

end 


NAC03300 

NaCO 33 1 0 
NftCO 3320 

NaCO 3.330 

NAC033A0 
NaCO 3350 

— NAC03360 

NaCO 3370 
NACO 3380 

N ACO 3390 

NaCO 3400 
NAC03A10 

NACO 3420 

NAC03430 
NACO 3440 

NaCO 3450 - 

NACO 3460 
NaCO 3470 

NAC0 3480 

NAC03490 
NACO 3500 

— NAC03510 

NACO 3520 
NaCO 3530 

— NaCO 3540 - 

NACO 3550 
NAC03560 

*AC03570 

NaCO 3580 


SUBROUTINE TAPER < IoEAO, TRITE-! 


C0hM 0N/8LK1 /Dl.lMl (37 3-1 ) 
COPMON/8LN2^DIJM2 (228^*) 
COPMCN/ULN3/OUK3(340T>- 
C0VMCN/«L K '* / 0HP4 ( 243-> ) 
C0pMCN/HLKS/nUM5 (402\ 

commcn/hlk 6 /oum 6 ( *o-i-i — 


C0^MCN/bLK7/DU87 ( I ) 

COf'McN/bLKR/OUPR (5) 

G-OmMcN / 8LK5 /Dt IR9 ( 6 ) 

1001 FOfJPA T ( 1 X i l OHTAPE R£*D .I5> 

1002 FORMAT ( lX» 1 lHTAPE wRfTE * 1 5 ) 

I F- ( 1 R E A 0 , E C . o )- GO-TO— 1 00 

REWIND 12 
DO 5(j J=1 » I W EAD 

— —WRI Tf (6* lOiil ) I READ — — 

5o RE aO { 1 2 ) ULM] t DUM2,DHM3,DUP*»0UM5»nUM6,DUM7,DUM8»DuM9 
RETURN 

— 1 00— I F- ( IR I It ,EC . 0 1 — ft E-T U R w — 

IF (lPITE.EC.-1 ) rEwImD 12 
WRI]E (6»1002j IRItE 

WRITE 1 1 1)._DUm 1. »DUM2 * rvUM3 *DUM4 •DUm5 ,DUm&-*DI>M7 *DUP8 »DUM9 

RETURN 


END 


-T 000 10 

T 00020 
T 0 0 0 3 0 
"T — 0 00*0 
T 00050 
7 00060 

-T — 00070 
T 0*0 8 0 
T 0*090 
T — 0*100 
T 00110 
T 00120 
-7 — 0*130 
T 001*0 
T 0*150 
-T— 00160 
T 00170 
T 00180 

-7 0*190 

T 00200 
T 00210 


T 0*230 
T 002*0 
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— supwpuT inE- point 

COV^L'TES ALL INTgRIOo AND 8OUN0APY POINTS 

COV^CN/tiL^^NC (B) •N’CfB) *MCl * NC2 *NCt . NC4 ♦ NC5 ♦ NC6 « NC7 ♦ NC8 * PC \ «MC2 

,.l_,MC3.MCA.MCP*Mt6,MC 7 ,Mc«.NKEG<8) ,No.C (2) »MMC (80»2) «nMA* ,MMx 

2 VGAMpA.GA.GB.GC.GO.GrtGF.X (40.81 ,Yf 19.8) «XXP( 130) 

3*YYP(l3d*i9> .Hm»XE. Yc.YAtXCfRn , Rn*r>MFLO» TT * CC .E* »Pl I 
-A.SIMhE (20 ) « COSThE (27) ,P <20 ♦ 19) .LSvMtLA.D* (8) .DY (8>- - 


C0f-MCN/BL'<2 / P ( I5n* 19» *U(150«l<5).Vfj50*19),S(150*l 9 )*PI (150*19) 
1»UI C 150 • 1 •* » *vi (150. In) ,SI (150*19) »mS(8.2) .NF(8.2) 

-2, MS 18.2) »MP 18.2) t TIMr«OT»K.J,XOTT(P) ♦ 3lNF-.QlNFN.KDlVS.DTS- 


C0MMCN/HL^3 / YET ( 1 3q ♦ T 9 ) .XCS<130) . XA » X 1 , X2 . X 3 . YO ♦ Y2 . ALP ( A > *DD (6 > 
l *8E 1 (* > 'LSLM.hU < 13() ) .Hl ( 130 > .hC < 1 3 ‘) ♦ HUPS (130) .NLPR U30) .HqPR ( 1 30 ) P 

— 0 In£nSI on -C OLfT(flf). 1 9’, ♦UOLO <80 *1 9 )-»vOL£>( 80 H 9 >-* SOLO (80.1 9 ) P 

NR1=NC1 

N«2=NREG (2) *NC2 

— NR-3 = NPEG ( 3 ) ♦ 1 
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105 . 


NR7=NREG (7) *1 
NL8= (NC8* 1 ) / 2 

-NM8=AREG (8 1 *NL8 : 

DO 2g L=l*2 
LRE6=3 

_LEL(.L .EQ . 2)-.A-_REG=.7_ 
NFlNsNC (L h EG) 
MFlNaMC <L r E°) 

;.00_16-N=1 t\KlU 

NM= (l-1) *N«- Ax*N 

nn=^reg <l r eg> *n 

.0 Q_1 6- M s l. IN 


1 0 
2 o 


POlO (NM.M) = p (NN. m) 

UOlO (Nm.M/ sU (NN.m) 
-YOLD(NM.Ml *V (NN. m). 
SOlO(Nm.m) =S(NN,M) 
COnTjNUE . 

.. 0 Q_ i 0 0 J_L M A c = L, 2 

00 4fi5 LREG = 1 ,7 
NSTA=NS<L H EO,LMAC> 
-NFINsNFILR^iLMAC). 
MSTAsMS (LRP-^tLMAc) 
MF I N r MF (L^E u .’*LPAc ’) 
-Dxx=ox (LREGl 


OYY = OY <lpeg ) 

IF (^SYM.NE. I ) GO Tc T 02 
IF (LrEG .G I . i . AND .LfiEft.Ll.*6l- 
C0\].INUE 

DO 400 M=MSTa,mFiN 
. L = H “ 1 * L K A C 


_GQ T_q 4 pS- 


IF(L.EQ.1)L=2 

IF (L.EQ.MC (LREG) *1 ) L4MC (LPEG) 

_YY = T iM.LREG) 

DO Ago N=NS f A , nF IN 
IF (LCEG.E^.P) ftO TC iOO 

.NNsNkEG 1 L r ^ ^ ) 

IF(Lsym.eQ.I) go TC Tor 

IF {Kn.EQ.NRI .ANO.M.E n.l) 


.IF.(NN.E(j.Nfil»ANO,M.En.l>- 


GO 

-GO- 


TO 4f)0 
.TO-AoO— 
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—oooio 
00020 
0 00 30 
— 000A0 
00050 
00060 
— 00070 
0O080 
00090 
'00100 
00110 
00120 
—00 130 
P 00140 
P 0O150 
P— 0O160 
P 0O170 
P 0 i) 1 8 0 
P 00190 


0O200 
0 02 1 0 
00220 
00230 
0 0 24 0 
— P — 00250 
P 00260 
P 00270 
— . P — 0 0 2 8 0 
P 0*290 
P 00300 
—P—00310 
P 00320 

0 0330 
-00340 
00350 
00360 
-00370 


P 

-P- 

P 

P 

-P- 

P 

P 

-P- 

P 

P 

P- 

P 

P 

-P_ 


P 

P 

P- 

P 

P 

-P- 

P 

P 

P- 


00380 
00390 
00400 
00410 
00420 
00430 
00440 
00450 
-00460 
00470 
00480 
00490 
00500 
00510 
0 0520 
0 0530 
0 0540 
00550 



IF(NV.EU.NR<;.AND.M.E«.mC2>-.-QO-JO-4Ao 

IF <*N .EO.NRJ.ANn.M.Eo»MC3) GO TO 4*0 
108 CONTINUE 

I=NN-1 *LMAC — 

IF (LyAC.EU.ii.aND.N.Eo^NC (L«EG) ) I=nN 
XXsMn.LPEG) 

GO— Io ( 1 1 0 . 1 20 »4 30 0 ♦4 50-4-160 *4-7 O-v-^LRFG — 

C coefficients 

110 csxn , = i. 

ety p =- 1 . : 

ETxP=0. 

GO TC SCO 

12 0 _CSxtsl» 

ETyPsI ./YA 
ETxPsO. 

GO— T 0-2 0 0 — 

1 3o CSxF=l./XC 

ETyR= W./hl (NN) 

E-TA = Y Y.p <NN,M> *ETYP 

ETxP=-ErA<»MPR (NN) «Er yp 
GO To 2 00 

IAO-CS xF XE-jcc ) 

£TyP=1 ./ (HL (NN) -HC (N,,) ) 
ETA = (YYP(NN»M)-MC(NN'| )*ETYP 

— — E-TxP-a META -.1.*..) oHCPfi < m N)«E-TA*hlPP-(Nm4-M»ET-YR 

GO TC 200 

1 So csxF=-i. 

- E.IyPs -U2.C YEthc ( NN J ) 

ETA=(YYP(NN«MJ-MC<NN, )*ETYP 
ETxP= (ETA-t . ) *hCPR (N m) *ETYP 

GOT 0—2 10 — 

160 CSxP=-l . 

ETyPs-1 . 

EIxP = KUPR CNN) 

GO TO 200 
17o CSXP=1./XE 

£.TyP=-L, 

ETxP=HUPR (NN) 

200 GO TO (225,250) ,lMAC 

C DERIVATIVES — 

225 PY=(P(NN,L)-P(NN,L-1\ J/OYY 
UYs (U (NNfU -U(NN,L-li >/OYY 

V-Y* l-V (NN. t L) -V (NN,L-I > ) /Q-¥¥ 

SY= (S (NN,L) -S (NN.L-l, ) /OYY 
PXs (P ( I *M) -P ( 1-1 ,M) ) /DXX 

JJX = -(U M-*M) .U(I-1^M)-) V 0XX 

yX= ( V ( I ,M> -V ( T-l ,M) ) /DXX 
SX= (S ( I *M) -S ( l-l ,M) ) /DXX 

P.P.s P (NN j M 1 

U U = U (NNi M) 

VVaV ( NN 1 M) 

SSa^ ( NN t M ) — — 

GO TC 300 

25o PY=.«PI (NN,L)-PriNN,L.l ) )/OYY 

UY* lUl-tNNil >--UI (nN.L^1M/DYY 

, 108 


P — 00560 

P 0*570 
P 015BO 

P 00590 

P 00600 
P 00610 

p 00620 

P 00630 
P 01640 

P 0O650 

P 00660 
P 00670 

P — 00680 

P 0O690 
P 0O700 

-P 00710 

P 0O720 
P 00730 

P-00740 

P 00750 
P 0o760 

P — 00770 

P 00780 
P 00790 

*>— 01800 

P 00810 
P 00820 

P 00830 

P 00840 
P 00850 

P — 0 0860 

P 0O870 
P 01B80 

P — 0 o 8 9 0 

P 01900 
P 0 09 1 0 
— — P — 01920 
P 0 0930 
P 00940 

P — 01950 

P 00960 
P 01970 
— — P — 01980 
P 01990 
P 01000 

P 0 1010 

P 01020 
P 01030 

P 01040 

P 01050 
P 01060 
— P — 01070 
P 01080 
P 01090 
P—01100 



VYo < V I ( NN » Ll-V I ( NN « L . 1 y ) /&Y Y 

SYb<SI <NN*L>-SI (NN«l. 1) l/OYY 
IF (LBEG.Nt.-J) GO TC ^55 ‘ 

IF (CMAC.NE.i!) 60 TO ^55 

IF (N.EO. I ) Go TO 255 

ACH =SORT((<J(nn .M)o«2*V(NN »M) **2j / (GAMMA«EXP (P (Nn tM)/GA.S(NN 

14/GAMMA) ) ) 

IF (Ach.GT. 1 *01) I=NN 
255 CONtlNuE . 

px= (p i-< i . m )- _ p i-(T— i .Myi/oxx 

UX* till (I tMl-UI (T-l »M> ) /oxx 
VX=<VI (I»M)-VI (T-l.M - , I/OXX 

SXs tSI ( I * M ) “ S l-M-r 1* M-V4-/0 X-X 

PPs^I (NN.M) 

UUsUl (NN.M) 

VV = V I •< NN.M) 

SS = ;?T (NN*M) 

P0=P(NN*M) 

UOsU(NN.M) 1 

V0=V(NN.M) 

SOrS(NN.M) 

30 C — CO K T I N U E 

XP=XXP(MN) 

YPsYYP(NN* w > 

A A s C c X P> *x C ' 1 NN )- 

AB=ETXP*YET (NN.M) 

AO=Et YP*YET (NN.M) 

AF-=Ou«AA__ : 

AG = C(j»Att* V V“AP 
AT=£XP (PP/GA*SS/GAMMy ) 

AH=0 . : 

L^X = 1 

It <ABS(YP> .LT.l.E-s>(_AX=0 

M«EG = ] - 

IF (LWEG.Lt.i .0R.l_RrG,GE.6) MREG=? ' 

IF(LA.NE.l) GO TO 3 1 ‘ 

IE. (Lax .Eq. o.)_ AH = AO« : . . .... 

IF (LaX.NE. ) AH = VV/Yo 
C EQUATIONS 

3 LO_P. T a - ( A F * P x. * a G * PY.tGA M u A.* (,AO?.UX.t A B *U.y> AD * VX. A H )4 . 

UT = - ( AF«UA*AG«ijY*AT* ) AA*PX*AB*PY) l-XOTT(LwAC) 

VT = * ( AF»vX*Ag*vy*AT<**D«py) 

£Tj»JLC AF_*SX ♦. Aq«s.Yj_ __ 

GO TC ( 350 « 375 ) »l_MAC 
35 t PI (Nk,M) =pp*pt«ot 

UI((NN«M). = UL*UT«DT : 

VI (Nn .M) =W»VT«OT 
SI (NN «M) =SS*ST*DT 

IE(MrEG.E 0-1 .ANn.M.E^l)_G.a_-_ZO_3A0 

IF (McEG.EU. i ,AN0.M,En.MC2) GO TO 3*2 
IF<MREG.EU.2.AND.M.E*.l) go TO 364 

GO .TO 4 JO — ! ! 

36n HPpIrHCPR (NN) 

GO 1C 370 

362 I F_l L R E G . E l±. 2) . _G 0 _1C_ a04) 


P — OHIO 
P 01120 


*M 


—P—0 11 30 
P 01140 
P 01150 

-P 01160 

P 01170 
P 01180 
— P—01190 
P 01200 
P 0 1210 
—P—01220 
P 01230 
P 01240 
— P — -01250 
P 01260 
P 01270 
— P — 01280 
P 01290 
P 01300 
~P — 01310 
P 01320 
P 01330 
-P— 01340 
P 01350 
P 01360 

-P 01370 

P 0)380 
P 01390 
~P — 01400 
P 01410 
P 01420 
_P.__014 30 
P 01440 
P 01450 
_P — 01460 
P 01470 
P 01480 
-P — 01490 
P 01500 
P 01510 
-P— 01520 
P 01530 
P 01540 
-P — 01550 
P 01560 
P 01570 
— P-- — 0 1 500 
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HPRlsHLPRlA.Nl ' 

GO <0 370 

364 IF (CrEG.EO.I.AMO.LSYm.NE.I) GO TO 400 

_ HPR, a .sHUPR < NNJ : 

37c SOn=SQRT(l .♦HPPI««2) 

T A 1 * I « /SQW 

TA2 = k.PRI»Tal 

VWaUU*TAl*W*TA2 
VWT = UT»TAl *V t*TA2 

V tet-I = V W ♦ V W TJ* U T- 

UI (N\«M1 =V*I*TAl 
VI (Nn ,M) =VwI*TA2 

: — GO— lo — * ! “0- 

» 375 P (nN,m) =, 5« (PP*Po*PT»OT) 

U <NN*M) = ,5« <UU*IJ0*UT*0T) 

V C N N ♦ M ) =.. 5 <» ( V V ♦-V O ♦ V T * 0 T-) 

S (NN»M) =.S« <SS*SO*ST*OT> 
IFcMreG.EO.C.AND.M.Ea.I) go TO 380 

IFCMREU.EQ.i.ANn.H.E/>«MC2V-60-T0-3«2— 

IF (MREG.EU.ii. AND. M.Ert.l) GO TO 384 
GO To 3S»5 

3flj HPrIsHCPR-LNN I 

GO IC 390 

382 IF (LREG.E0.2 ) Go TO 395 

: HP.Ri.sHL PR.IA.N4 

GO TO 390 

384 IF(CpEG.EU.1.AND.LSy».NE.1) GO TO 400 

HPR AsHUPR (NN) __ 

390 SQR=cQHT (1 .*HPRI**2) 

T A l a i ./SrjR 

TA2 = hp«I*TA1 

VW=Uu*TAl *VV* ta2 
. VWO = UO*TA1+Vo*tA2 

V W-T = u T .* T A l .+.V.T-# T Aj> 

VWN = ,5* (VW*Vwo*VWT«Ot) 

U (NN«M) =VW\*TA1 

V-INN.MLsVWN*TA2— 

395 CONTINUE 
400 CONtlNUE 

4.05— CO f. t INUE— : 

C CUS^ point 

IF (Lsym.EO. 1 . ANO.LHAr.EO.il GO To 500 

IF_(LS YM .eO«I » ANO»LMA(L«EO»2-1 — GO— TO— p 00- 

NA1=NCI-1 . 

MA3=PC3-1 

PK8 =7 ; 

YP= I YP <NRl 1 1 ) 

AHsO , 

G 0 _. Tc — (4 1 0 » 420 ) .LMAC- , 

4 1 0 OXP s XXP<NRt)-Ro*R(NLn»MM8) 

0YP = yTP (NR3 *MCD -YYP',NQ3,mA3) 

RXP= (P (NRl . 1 ) -P ( nM8 » mMA ) ) /CXP 

UXP= (U (NRl « 1 ) *U (NH8,,.M8) 1 /CXP 
VXP 3 (V (NRl » 1 ) ♦ V ( nM 8 « uM() ) ) /CXP 

$XP=(S<NRl»l)-S<NM8t.,M8> )/DXP- 
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— P 01590 

P 01600 
P 01610 

— P — 01620 
P 01630 
P 01640 

— P - 01650 
P 01660 
P 01670 
— P — 01600 
P 01690 
P 01700 
-P—01710 
P 01720 
P 01730 
— P — 01740 
P 01750 
P 01760 
—P—01770 
P 01780 
P 01790 
— P — 01800 
P 01810 
P 01820 
— P — -01 830 
P 01840 
P 01850 
— P — 01860 
P 01870 
P 01880 

P 01890 

P 01900 
P 01910 

— P 01920 

P 01930 
P 01940 
—P—01950 
P 01960 
P 01970 

— P—01980 
P 0 l 990 
P 02000 

— P 02010 

P 02020 
P 02030 
— P — 02040 
P 02050 
P 02060 

— P—02070 
P 02080 
P 02090 

— P -02100 
P 02110 
P 02120 
— P— 02130 





pyps (P (NRl ♦ l ) -p (nR 3 t«.A3) » /OYP 

UYP= (UINRl. 1 >-|)(NR3 »uA 3) )/0YP 
VYP a (V (NRl » 1 ) -v (NR3»«/A3 j ) /DYP 

SYP a <S (NRl « 1 )-S (NRJ»..A3> )/OYP 

PPaP (NR 1 » 1 ) 

UUsU (NR 1 « 1 ) 

VVsV (NRl. 1 )- 

SS=5(NRltl) 

GO Tc A3o 

-4.2fl-~0XP = XXP (NRl 1 -XxP-(NA-l', 

DYP = yYP(NR 1 *2) _YYP (Nol . 1 ) 

P Xp = ( P I (NRl * 1 ) « P I < N A 7 . 1 ) l/OXP 

UXp s (UI <NHl-*-lY-UI <NAr» l-H /DXP — 

VXP* (VI (NRl * 1 ) -VI (NA?#1 n /Dxp 
SXps(SI(NRl»l)-Sl<NAT.l))/l)XP 

PYP=(PI <N«l *2)-- PI (NRTh H/D-YP 

UYP = (UI (NRl *2) -UI (NRT » 1) ) /OYP 
VYP= (VI (NRl *2) -VI (\R7» l J) /OYP 

SYp *-( SK nRI * 2 l-SMNRf* 1-H-/DYP 

PP=Pt (NRl « 1 > 

UUaUr (NRl * 1 ) 

VWsV.I (NRl-U ) 

SS=?I(NR1*1> 

POoP (NRl . 1 1 

U0sU4NRU.il 

VOsV(NRUI) 

S0 = i>(NRUI) 

ASH-.! E.(L A » Nt. ll_GQ.-TO. AA-g 

AHsVv/YP 

44c ATsJrYP CPP/GA*SS/GARM; ) 

PJT=- ( U U ? R AEA V.V * P Y P>g ; MMA«4 U XP*.V„yp.* aHU 

UT=- (UU*UXF*VV*UYP*At«PXP«XOTT (LMA r) ) 
VT=- (UU*VXP*Vv*VYP*At*PYP) 

S Is.- ( UU*SX.= ♦VV«SYP) 

GO To (450 *469 ) *LMflC 
45o PI ( fY K 1 * 1 ) =pP + pT *oT 

— : — U I.(Nc U Lj sLU*in:«oT 

VI (NRl . 1 ) =VV*vT«oT 
SI (Nr1*1)=SS*sT*dT 

EC C-P LCNRU.l 1 

UCC = UI (NRl « 1 ) 

VCC = Vl (NRl *1 1 

S C C .= SLi NR l.t.1 1 

PI ( NR7 * 1 ) =PCC 
UI (Nr 7* 1 ) =l-Cc 

.V.I<Nfi7.1) =7CC; 1 

SI ( N R 7 * 1 ) =«0c 
PI (Nc2»MC^l = PCC 

UI (Nr2..MC2> =UCC 

VI (NR2«MC2) =VCC 
SI (Nr2* MC2) =SCC 

P-I (Nfi3tMC3l =PCC 

UI ( Nc 3 t MC3 ) =UCC 
VI (Nr3*MC3) =VCC 

SI (NR3.MC3) =SCC— 


P 021 40 

P 02150 
P 02160 

p 02170 

P 02180 
P 02190 

P - 02200 

P 02210 
P 02220 

P- 02230 

P 02240 
P 02250 

P — 02260 

P 02270 
P 02280 
— -P — 02290 
P 02300 
P 02310 

P— 02320 

P 02330 
P 02340 

P—02350 

P 02360 
P 02370 

P — 02380 

P 02390 
P 02400 

P 0 2410 

P 02420 
P 02430 

P—02440 

P 02450 
P 02460 
p_. 024 70 

P 0248(1 
P 02490 

P — 02500 

P 02510 
P 02520 

P—02530 

P 02540 
P 02550 

P — 02560 

P 02570 
P 02580 

P — 02590 

P 02600 
P 02610 
— P — 02620 
P 02630 
P 02640 
— P — 02650 
P 02660 
P 02670 
— P—02680 
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GQ_.l0 -SO 0 

460 P (N?j * 1 )=.«=* (PP*p04Pt« 0T> 
U ( fvRl « l ) = .5* (1 )ii*uO*Ut*OT) 

V { f*R 1 « 1 ) s • 5 * (W*VO*Vt®DT) 

S (N«l ♦ 1 ) «.5* (5S*S0 *St*0T) 
PCC S P (NR 1 * 1 ) 

UCC=U (NR1 1 1.) 

VCC 3 V (NR1 * l ) 

SCC = S (NR1 * 1 ) 

PI N R 7. *_1 ). = PC C 

U(kR7*1)=UCC 

V (nR 7 . 1 ) =vcc 

S.lKRlt L).= SCC 

P (N'R2«MC2) =PCC 
U (nR2»mC 2) aUcc 

_.V.INR2j MC2 ) 5 V CC 

S (NR2.MC2) =SCC 
P (nR3«mC3) =PCC 

LMMR3<mC3) =UCC — : 

V (NR3*MC3) =VCC 
S (nR 3»MC3) a See 

GO -T-C— 8f#-0 

C FRE£ STREAM 

500 CONTINUE 

00-510 -Mai,. MCI 

51c UI (1,M)=QINF 
DO 520 Nsl.Nci 

520 4J I (N,MCl)=ClNF 

NN2=N p EG (2) *1 
00 530 M=ltMC2 

530 .UI (NN 2 .M) sClNF- 

NN5=NREG (5) ♦NCS 
00 ?40 M=1iMc«5 

5<tO_U I. CNN 5 « M ). =G 1 NF 

NN6=NREG (6) *NC6 
00 550 M= 1 • M C6 

55.0 — U I.(NN6 » M>-sClNF 

00 ?60 N=l,NCft 
NNsNREG (6) ♦ N 

5 6n UI (Nn , MC6 ) = V.I NF. 

00 ?7q N=1«NC7 
NN = "REG <71 *N 

— — 57 0- U I ( *N » MC7 ) =U I NF 

600 CONTINUE 

C INTERFACE A 

I F.(L$Y v . E Q . J.l _G 0— T C-a 1 5 

NCiL=NCl-l 
00 610 N=l,Nci(_ 

NNsNrEG.(2 )-♦ N 

PI (N» 1) sPlINN, RC2) 

UI (N, 1 ) =Ul <N n,mC2) 

' VI (N,l)rvI.(.NN f KC2l : 

6 1 C SI <N« 1>*SI (NN»HC2> 

C INTERFACE b 

615-CONT INUE — 


9 — 02690 

P 02700 
P 02710 

— 9 02720 

P 02730 
P 02740 

— P 02750 

P 02760 
P 02770 

P — 02780 

P 02790 
P 02800 
-P—02810 
P 02820 
P 02830 

r-P — 02840 

P 02850 
P 02860 

P — 02870 

P 02880 
P 02890 

P—02900 

P 02910 
P 02920 

P—02930 

P 02940 
P 02950 

P — 02960 

P 02970 
P 02900 

P -02990 

P 03000 
P 03010 

P — 03020 

P 03030 
P 03040 

P — 03050 

P 03060 
P 03070 

P — -0 3080 

P 03090 
P 03100 

P — 0 3 1 1 0 

P 03120 
P 03130 

P — 0 3140 

P 03150 
P 03160 

-P—0 3170 

P 03180 
P 03190 

P — 0 3200 

P 03210 
P 03220 
; P — 03230 



! NN 2 SN R E G < 2 > > NC 2 — — 

NN 3 = NREG { 3 ) ♦ 1 
00 620 M=liMC2 

Pt<NN3»M)=cX <NN2,M> 

UI (Nn 3*K) si, I (NN2,M> 

VI (Nk3'.M-)*vI (NN2,M) 

62«-SI <*N3.M)=SX ( NN2 *M I 

C INTERFACE C 

NN 4 =NREG (4) *NC4 

1 F < >r S UP . G T . 0 | -g O-T C -*3?- 

00 030 M=1.*C4 
PI ( NN4 » M ) s-2,3 

&30 — SI ( Nf,4 » M ) 

‘ GO TC b38 

632 NN4 L i =NN4 “ l 

-NN4L2=NN4-2 

oo 635 m=i,mc 4 

P I ( NN4 * M ) =2 • *P I (NN4|_T tM)-Pl<NN4|_2«M) 

-UI (Nn4»M) s2« #U-I (NN4L-v»M) -UI-<NN4L.?*m+ 

VI (Nn4 ,M) -2 «*V 1 (NN4LT *M) -vl (NN4L2 »“) 
635 SI(.NN4 »M)s 2.»SI (NN4LT »M) -Si (NN4|_2 *m) 

63S--C0nT IMUE ; 

C INTERFACE E 

NN 6 =NREG ( 6) *1 

-HtTf .= VR E G i U ♦ N C-7 : 

DO 650 M=1*MC6 
PI <NN 6 *K) =FI (NN7 # M) 

u J.( Nk "6 • Ml -U U NN7-*M > — 

VI (NNft«M)=Vl <NN7,M) 

65C SI (NN 6 «M)s«I (NN7»M) 

_£ -INTERFACE F 

NN7=NPEG(7) tl 
DO 660 M=1 ,MC7 

PJ INN7.I.M) =EUnc1.M) 

UI (NN7,M)slI (NCI ,M) 

VI (NN7»M) =Vl (NCI «M) 

660—S I <Nn 7.M) ssl.<NC-l iM) 

C INTERFACE g 

NN3=NR£G (3) *NC3 

MN.4 = NREG.(A)AJ 

00 670 M s 1 » M C 3 

PI ( Nn 4 * M ) = R I ( NN 3 * M ) 

U I.(N\4 » M).=l I (NN3,M)_ 

VI ( N N 4 * M ) =VI (NN3,M) 

67c SI (Nn 4 »to ) *S I (NN3*M) 

GO -T C-_-l 0 0 0 : 

C FREt- STREAM 

BOO CONTINUE 

00-81 0 -M=1» MC-1 

Bin U(lfM) sOINFN- 

00 820 N=l,Nci 

82n__U (MP.C1 ) sUIUFN 

NN2=NREG (2) *1 
00 830 M= 1 .MC 2 

; 83c_U.CNN2 iMI.sUINfn 1, 
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-P — 03240 
P 03250 
P 03260 
-P — 03270 
P 03280 
P 03290 
-P—03300 
P 03310 
P 03320 
— P — 03330 
P 03340 
P 03350 
—P — 0 3360 
P 03370 
P 03380 
— P-03390 
P 03400 
P 03410 
-P — 03420 
P 03430 
P 03440 

-P (13450 

P 03460 
P 03470 
-P— 03480 
P 03490 
P 03500 
— P — 03510 
P 03520 
P 03530 
-P— 03540- 
' P 03550 
P 03560 

P 03570 

P 03580 
P 03590 
-P — 03600 
P 03610 
P 03620 
-P — 03630- 
P 03640 
P 03650 
-P — 03660 
P 03670 
P 03680 
-P — 03690 
P 03700 
P 03710 
-P— 03720 
P 03730 
P 03740 
-P — 03750 
P 03760 
P 01770 
-P—03780 



NN5 = NREG.(5)j»NC5 

00 840 M=1 »mcs 

8*0 U (NN5 *M ) stJINFM 

NN6*NREG (6) *NC6 

DO 850 M=1*MC6 
850 UlNNfttM) sUINEN 

OQ-860 -N=l »NC6 

NNsNpEG lb) * N 
860 U (NN»MC6) aUlwFN 

00—87 '>-N=l »NC7 

NNsNREG < 7> *N 
870 U(NN.MC7) =UINFN 

— aoo— continue ; — 

c interface a 

IF(f-SYM,EO.i) GO TC ol5 

NClL-sNCl “ I 

DO 81 0 N=1*NC1L 
NNsNREG (2) *N 

P < Nn « Mt 2 >-=P- < N r H 

U(Nn»MC2J=U(N*1) 

V (Nn »MC 2) =v (N, 1 ) 

9lo — S (N\,MC21 = !i (Nil) i 

c Interface « 

915 continue 

— NN3 s NREG121-*NC2 

NN3=NREG<3) *1 
DO 920 M=1»MC2 

E_(NN2iM.).bP_(Nn3,Mj. 

U(NN2»M)sU(NN3»M) 

V (NN2»M) sV (NN3*M) 

920-S (NN2»Ml =S INN3 *m> 

C INTERFACE C 

NN4=NREG<4) *NC4 

IE.<LSUR-.GT-.U > -60_TC-^3a- 

00 930 Msl*flC4 
P ( NN4 t M ) s-2 • 3 

93a-S (NN4 * M 1=0 — ; : 

GO TC 938 
932 NN4Ll=NN4-l 

N N4L2 =NN4^ 2- 

DO 935 N=ltMC4 

P (NN4»M)=2.*P (nN4lT»M)-P (NN4L2*m) 

u_- <NN4*M) =2.*u -(NN4L-T-* MT^-U -1NN4 L?*m> 

V (NN4iM)=2.*V (NN4 l7»M)-V <NN4L2*m) 
935 S (NN4 *M) =2.*S (NN 4LI n<) -S (NN4L2*M) 

SI 8_C 0 N T I NUE... 

C INTERFACE E 

NN6=NREG(6) ♦! 

NN7=nREG (7-!*NC7 

00 9so Mal,MC6 
P (NN7*M) = P ( NN6 * M ) 

UXNN7 »H1 =0 (NN6iM) 

V (NN7*M) = V (NN6iM) 

95c S <N'N7»M) sS (Nng.M) 

-C INTERFACE -E : 


— P — 03790 
P 03800 
P 03810 
— -P — 0 3820 
P 03830 
P 03840 

P — 03850 

P 03860 
P 03870 
— P — 03880 
P 03890 
P 03900 

— P— 03910 
P 03920 
P 03930 

P — 03940 

P 03950 
P 03960 
— P-03970 
P 03980 
P 03990 

— P — 04000- 
P 04010 
P 04020 

— P — 04030 
P 04040 
P 04050 
— P — 04060 
P 04070 
P 04080 
— P — 04090 
P 04]00 
P 04110 

— P 04120 

P 04130 
P 04140 
— P — 04150 
P 04160 
P 04170 
— P — 04180 
P 04190 
P 04200 
—P—04210 
P 04220 
P 04230 
— P— -0424 0 
P 04250 
P 04260 
—P—0 4270 
P 04280 
P 04290 
— P — 04300 
P 04310 
P 04320 
—P— 04 330 
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NN7=NREG(7)-M 

00 ;60 M=1«MC7 
P(NCi,M)=P(NN7,M) 

U ( f,C i *m ) =U (NN7 — 

V (NCI ,M) = V (NN7»M) 

960 S (kCi »M) =5 (NN7»M) 

— C— INTERFACE- G — 

NN3=NREG <3) *NC3 
NN4 = N REG ( 4 ) ♦ 1 

00 -^70 m=1-,MC3 

P(NN3iM)=P(NN4tM) 

U (NN3 »M ) =U (NN4 «M) 

V (NN3»M) »V (NN4 t M) — 

97c S (NN3 iM) =S (NN 4 f M) 

1000 CONTINUE 

00 1 o 2 u - L= 1 * 2 

LRE0=3 

IF ( L • EQ • 2 ) LREG = 7 
NF-INoNC-ILREO.) 

MFlN=MC (LR£^) 

00 i o 1 o N=1 *NF IN 

N M = .( L ^-1_ ) * NM A X * N____ 

NNrNREG (LREG) *N 
00 1 fl 1 n. M s 1 * MF JN 

Pl.(NN.tM)=PCLD(.NM,M) 

UI (Nn ,M) =UCLD (NM,M) 

VI (Nn*M) =VCLD(NM,M) 

l-O-l-O — S I t f^N * MJ sSC.LD.f NM* M ) 

1 020 CON ' INUE 

00 Io3G L=l*2 

LRE G = 3 ________ 

IF<L.EQ.2) URFG=4 
MFINaMC (LREG) 

NFlN = NC.(LREOl.,i 

00 1 028 M=1 *MFIN 
00 lfi24 Ns 1 * NF I N 

NN = NREG (|.REG>-*N. ....... 

ACH =S0RT((O(nN »M)*i2«V(NN *M) **2) / (GAMMfl«EXP (P (N n «M)/GA*S(NN iH' 
1 )/GAwMA) ) ) 

I F. { A c H . L II. 1 j*. 0-U—GQ— In— 1 0 24 

N1=N*1 

NL1=NN*1 

ACH I = SORT 1(W(NL1*MU»o24\/-(NLL« M G A MM a « E-X P( P- ( N|_ I * M ) /G A ♦ S < NL 1 ♦ M - 

1 ) /GAMMA) ) ) 

IF(A C H1 .LT.I.oI) Gc t 0 1026 

GO -T.C - 10 24 

1 626 N2=**2 

NL2=NN*2 

P. ( N h 1 • M ). s t FJ NN , M ) ♦ F.CiiL2* M ) 1/2 . 

U(NLl ,M) = (L <NN.M) *u ( mL2,M) ) /2, 

V(NL1»M)=(V(NN,M)*V UiL2»M) )/2. 

S ( N L 1 » M L=_(JJ.NN ( (jL2_» M ).) _Z2^__ 

K-24 CONTINUE 
1 628 CONTINUE 
LQ3L£L_.C0N?.INUL._ 


-P-04340 
P 04350 
P 04360 
P — 04370 
P 04380 
P 04390 
— P — 04400 
P 04410 
P 04420 
— p- -04430 
P 04440 
P 04450 
— P -04460 
P 04470 
P 04480 
— P—04490 
P 04500 
P 04510 
— P— 04520 
P 04530 
P 04540 
— P - -04550 
P 04560 
P 04570 
— R — 04580 
P 04590 
P 04600 
—P—04610 
P 04620 


T15 



_NN3=NREG.(3UNC3. 


NN4=NRFG(4U\ 

DO 1 0 32 M* 1 .* MC 3 
P(nN9.M1sP(NN4.mi 

U (NN 3 »M)sU(Nna,M) 


V(NN3,M)sV(NN4.M) 


Lft32_S (KN-»,m) sS (NN4.M) 


. RETURN 

P 04630 

ENO 

P 04640 


___s__JSUBH.CUJ.lNt NOSE - — — 

COf'MrN/BLKI /NC (fl ) *NC1 » NC2 . NC "> ,NC4 .NC5 » NC6 « NC7 . NC8 . MC 1 ,MC2 

1 .MC3.MC4.MCP»ML6,MC7,Mca.NMEG<8 ) .NmC (2) «Mmc ( a 6 • 2 ) ♦NMAX.MpAX 

2.* GAMMA.t.GA.»OU..»GC «GO.»Gc.*GF »X (4 0 ,8 u Y t 1 9 .81 * *XP< 1 30 ) 

3 * Y YP ( 13ij» l*’* tHmXEtYn Y A , xC » Rl . RD . eMFLO. TT . CC . EP « P I 1 
4»SI ( UHE(20) *C0STHE(2*) ,«(20*19) »LSvM*LA*D)f (8) »0Y(8) 

C0mMcN/bLK2/P..< 150*JR'y • UO50.19UV( 150.19US ( 150.19+.PI (150.19>— 

1*UI ( 1 5 r * 1 V ) «vi (150*1(1) ».SI (150*19) »vS(8»2) .NF(8.2> 
2»MS<e«2)«MF(8,2),TlMr*OT,K»J*XOTT(5) tQTNF .(5INFN.KDlVS.OTS 
C.0 Y M C NZ9 L K 4 / E P 1 { 20. Ur P2 (2r.Jj£ P 3 .( 1 0 k ) . f.EP_4_(tO 05- . E P5. ( 3 8 ) . E P 6 1 3 a) __ 

1, NN a (2C) .NN2(?0) «M1 (^0) ,M2(20) «M-?C»8) »M4(38) .LI (100) *L3(100) 

2. H (lfUj . 12 UoOl.OP In8, 19). 8 (20) ,8 dR(2m 

OY Y 3 0 Y ( 8) :— 

DXX=CX (8) 

DO 5 o 0 KS=UKDIVS 

D.O-f (VO-.LMACAl., 2- — 

XACS:= (FLOAT < KS-2*LPAr> /FLOAT (KOIvSi ) « ( XOTT (2UX0TT ( 1 ) ) *X0TT < 1 ) 

. N S T A s N S ( H . L M A C ’) 

NFINsNF (8.LMAC+ 

MSTArMS (8tL M AC) 

MFINsPF (8.L M AC) 

00-i50~M=M5 TA.mFIN — : 

L = M"1+L , 'AC 
IF (M.EO.l ) L=2 

DO *5ft ..Ns NS J A « NF-I-N : — 

NNsiYREG (8) *n 
I=N h-l .LMAC 

IF-(4,EQ.NH£U.(8) .* 1 ) — = 1 — 

RRs'« (N.m) ' 

AAsi./(CC-e <N) ) 

AB = i . /PI-I— — 

AO= ( Y (Mi8) - 1 . ) «RPR (N) *AA 
GO Tc (25*^1 »L M AC 

25_P Y s Ip.LNN « L I. -P.t NN » L-l'v .1 /QYY 1 

UYs <U (NN.U “U (NN.L-l i )/OYY 
VYs ( V (NN.L) -V (NN.L-l » ) /OYY 

SYs lS,(NN.LJ -”S(NN.L*lV) /nYY 

PX= (P ( I *M)**P( 1-1 »M) ) /OXX 
UX=<U(I»M>- u (T-l,M) ) V DXX 

VXstv (I.M)-V (I-l.M) J /DXX 

SXs (S ( I*M)-S>( I-J ,M) ) /DXX 
PPcP ( NN * M ) ■ . 

UUsU (NN«M) - : I 

VVsV(NN.M) 

SSsS (NN.M) ' ° 


— NOSOOOlO 
NOS0T020 
NOS00030 
— NOS0O04O 
NOS0T050 
NOS0O060 
— NOS00070 

noSoooso 

NOS00090 
_J40S00100- 
NOSOOllO 
NOS00120 
-NOS0O130 
NOS00140 
NOS00150 
— NOSOO-V60- 
NOS0D170 
NOS0T180 
— NOS00190 
NOS0^200 
NOS00210 
--NOS0O220 
NOSO'1230 
NOS0O240 
-NOS00250 
NOS00260 
NOSOT270 
— NOSOOPSO 
NOS0T290 
NOS09300 
- NOSO-1310 
NOS00320 
NOS00330 
— NOS 0^340 
NOS00350 
NDS0T360 
— NOS0^370 
NOS0T380 
NOS0<i390 
— NOS0O400 
NOS0O410 
NOS0J420 
— NOS00430 
NOS00440 

NOSOiY4Sn 



GO-TC-ltO 

50 PYs IplINNiLl -PI (NN«L-1) )/DYY 
LlYs<UI <NN*L * “ 1*1 (NN» L»1 ) )/DVY 

' V_Ys.( V I (NN»L ) -V_I CNN«L.1UWD-YY_ 

S Ys (SI ( NN t L * *S I ( NN » L « 1 ) l/DYY 

PX=(PIU,M)-PI(I-1,m, )/0XX 

-UX= (Ll (I.M)-UI M-l,M'v)/n.XX 

VX= <vl ( I.M)-VI ( 1-1 ,M, )/r»xx 
SX3 ( s I (I.M)-SI U-1«M| J/OXX 
PE=E4-tNM*aJ : 


-NOS 0 O 460 
NoS 00470 
NOSOO 40 O 
-NOSO'1490 
NOS 01500 
NOSO 05 1 0 

^OSOo 52 o 
NOSO 0530 
NOS00540 
NOS00550 


UU=L/J <rgN f MT — — 

VVsVl (NN.M) 

SS = $I (NN.M) 

l€0-AE = Uf*4A*VV/RP*A0 

AF=Vy/PP«A8 

AT = i:xP (PR/G^SS/GAMMa ) 

AKs <C : U*5 I NT Mg (N, ) *Vv*^asTHF fN) 

PTs-(AE*PY+AF«P¥*GAMuA* < AA«UY+UU/Ro*AB*VX/RR*AD*VY/RR) ) 
I p < u A , E (J , 1 ) PT=PT-GAM„A*AK 

UTts- (AtouY*AF<njX_VV«»2/PP»Ar«AfropY v XACC*C0STHEM^ H 

VT = *” ( At 0 VY*Ar*VX + (JU*</V/OR*AT # (Aa*Pv*AQ*PY)/RR 
l-XACC*SI,MTHt (N) V 

S-T* M -( AE»SY-* a F+SX > 1 — 

GO le d 10.120) *I_MAC 
110 pi ( inn ,m) =pp*pt*ots 

: UI-(N\ . M ) =UU ♦UT«OTS ... 

VI (Nn.m) =VV*VT«r)TS 
SI (i.mm)=SS*ST«DTS 

— IF-(M # -NE , l-)-dO— TO— 1 50 

SQR=SORT ( 1 . ♦ (PPR (N) /o <N) ) *«2) 

TA2=1 ./SqR 

iai.=ta 2 *hpp. in.) /n (N» 

VW = IJu*TAl ♦VY*TA2 
VWT=UT*TAl * VT«TA2 

V-W N = V W *V W To U T s 

UI (NN »M) aVV.N*TAl 
VI <NmM)=VV.N*TA2 

GO -To 1 50- — 

12'! P <nN*M) =.5« (PP*P(NNi..) ♦PT*0TS) 

S <NN,M) =.S" <SS*S (NN ».,) ♦ST^DTS) 

IF-f.M.EQ. 1 ) _.GO— TO-U0 

U (NN.M) = ,5« <UU*U(NN ...» *UT*DTS) 

V (NN,M) = .5« ( VV*V (NN «..) ♦VT*0TS) 

GO--IC--J50 

14" SQR=SO«T ( 1 . ♦ (RPR (N) /a <N) )**2) 

T A2 = l ,/SqR 

__TAl=TA2*riPF-(N)/0 (N)_.: 

VW = U (NN.M) * T A 1 *V (NN,..) *TA2 
VWI=lU <, Ta 1 *Vv*tA 2 

VWT = L'T*TAl ♦ Vt»TA2 

VWN = ,5* (V*>*VWI ♦VWT«OtS) 

U(NN,m 1=V*N*TA1 

V (NN.M) =VW.'»*-TA2 __ 

15o CONllNUE ^^7 

200 CONllNUE 


-NOS00560 
NOS 0 O 570 
NOSO^SSO 
-NOS 0 n 590 
N0500600 
NOS00610 
-NOS 00620 
NOS00630 
NOS0064O 
-NOSOofeSO 
NOS00660 
NOS0O670 
-NOSOofeBO 
NOS 0 O 690 
NOS0T70O 
-N0S0O710 
NOSO 1720 
NOS00730 
— NoSO •'74 0 
NOS00750 
NOS0t760 
-NOSO->770 
NOSO 0 780 
NOS 0 O 790 
-NOSOoftOO 
NOS01810 
NOSO 1820 
N>OSOo830 
NOSO 0840 
NOS00850 
--NOS0 1860 
NOS0T870 
NOSO^BBO 
-NOS00890 
NOSO fl90 0 
NoSOoplO 
-NOS0O920 
NOS00930 
NOS0OR40 
-NOSO-iRSO 
NOS00960 
NOS0O970 
-NOS0T980 
NOSO °990 
NOSO 1000 



jC 


-INTERPOL A T UN- A-I- 0 UTcR-RO UNO A R Y R.C 

OEl s FLOAT <k 5 ) /FLOAT <*OlVS) 

MA8=PC8-1 

-KA8 a NCd-l 

DO 3 oO n= 2 »nah 
LL l a NN 1 ( N > 

U.2=NN2 (NL- 

JlaPl(N) 

J2=M2(,M) 

EPS 1 aEPl-INi : 


NoSotoio 
NOS 01020 
NOSOI 030 
NOSOIOAO 
NoSO 1050 
NOSOI 060 
-NOSO 1 0 70‘ 
NoSO l 080 

NOS01090 
NOSOI 100 


EPS 2 aEP 2 (N) 

NN=NREG( 0 >* n 

PlcOEL*(P(LL.l.jl )-PI ,LL 1 . jl) ) ♦PI (L( 1 • J 1 > 

Ul sUe L«_ tu LLUl^ j u -L LL !-♦ JUUU U Lf -1 • J 1 ) 

VlsOEL* < V <LL 1 , Jl ) -VI fLLl » Jl) ) *VI (Lj 1 *J 1 ) 
Sl=pEL*(S(LLl,jl ) -ST *LU . J 1 ) ) *S I (Li 1 * J U 

_P 2 * 0 EL«-(P-<El. 2 ^u?E-P I VLt? vj 24 -> *P HLE 2 » J?l 

U 2 = 0 eL* (U (LL 2 . J 2 ) -L I ',LL 2 « J 2 ) > *UI (Li 2 t J 2 > 
V 2 =UEL*(V(LL 2 .j 2 )-vr,LL 2 . J 2 ) >*VI<Li 2 »J 2 ) 

S 2 = 0 EL*- 1 S 1 LL 2 -* J 2 ) -S I <LL 2 ♦ J 244 ~*-SI-<Le 2 *J 2 -) 

PB=? 1 ♦EPSl* <P?-P 1 > 

UB=Ul ♦EPSl° (U 2 -U 1 ) 

UBaVl ♦EP.Sl* (V 2 -V 14 

S 0 = 9 \ ♦epsI* (S 2 -S 1 ) 

UC = y (NNfMAS ) *Cf)STHE (»|) -V{MNtMA 8 ) oStNTHE(N) 
VC=G 4 NNj>tAfl.).*SlNTHE-lMl*U(NNfMA 8 )-*CoS^tME- 4 N>- 

P (NNfMCrt) = c <NN,m* 8 ) ♦rPSP* (PB-P INN.-A 8 ) ) 

S (\N»MCb) *S (NNtMAS) ♦cPs 2 « ( S 8 -S { NNt wA« ) ) 

UAsUC*EPS 2 « tUB-UC) 

VA = VOEPS 2 « ( VR-VCl 

U (nN.mCH) =IA*C 0 STHE (k.) ♦VA«SlNTHE <N > 

- 3 jaO-J/-tKN»MC 8 X = -UA*clNTBEVN)-^UA*COSTHE-tM) — 

DO Ago L=l *2 
IF(L,EQ. 2 > GO tO 410 

LRE^=7 * ^ : 

NA A* 1 
NCC S 2 

GO— IC-A 28 

410 LREG =3 
NA A *NC 8 

nccsnC 8 =.i 

420 NNB=NREG (LPtG) *NNC (L\ 

NNA S NPEG ( 8 ) ♦NAA 

J<LC?KREG 18 1 ♦NCC 

DO 4 cq M= 1 ,MC 8 
MRb.Il -1 ) «M*'AX»M 

MM 3 S P 3 (MR 1 -^ i , : : 

MM 4 *P 4 ( HR ) 

EPS 5 =EP 5 (MR) 

:_EPS 6 = EP 6 (MR )__ 

P 3 = 0 EL« (P (NNB,MM 3 ) -Pt (NNB,PM 3 ) ) >pl ,NNB,MW») 
U 3 =DeL"(U(’>Nb«MM 3 )-u; (NNB*MM 3 ) ) fill / NN 8 « MM 3 ) 

V 3 »L>EL* ( V (NNb^MM 3 ) -V i ( NNO« MM 3 ) 1 * V I (NNB »MM 3 ) 

S 3 =pEL*(S(N N B»MM 3 )-S> (NNB.MM 3 ) ) ♦Rl(NNB*MM 3 ) 
P 4 * 0 EL* (P(NNB,MM 4 )-Pt (NMB.MM 4 ) ) *Pl <NNB,MM 4 ) 


-NOSOlllO 
NOSOI 120 
NOSO 1130 
-NOSO 1140 
NOSOI 150 
NOSOI 160 
-NOS01170 
Nosoueo 

NOSOI 190 
-NoSO 1 - 20 G- 

NOSO 1210 
NOSOI 220 
-NOS01230 
NOS01240 
NOSOI 250 
-NoSO 1260- 
NOS01270 
NOSO 1 280 
-NOSO i 290 
NoSO 1300 
NOSO 1310 
-NOS01320 
NOSO 1330 
NoSO 1340 
-NOSO l 350 
NOS01360 
NoSO 1370 
NOSO 1380 
NOS01390 
NoSO 1400 
NOSO 1410 
NOS01420 
NOSO 1430 
-NOSO 1440- 
NoSO 1 450 
NOS01460 
-NOS01470 
NoSO 1 480 
NOSO 1490 
NOSO 1500 
NoSO 1510 
NOSO 1520 
-NoSO 1530 
NOS01540 
NQS01550 
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_D4 = DEL* ( U {A'^B » MN 4 ) r.U.t-(NNB « KM4.Ll_»lll.(.NN3 »MMa) — 
VAspELMV <NNB,^M4)-Vi (NNR.MM4) > ♦ Vi ( NNH . MM4 ) 
S4=Uf L* <S <NNR,MM4 ) -Si (NNR ,MM41 ) *Sl (NNH.MM4) 

J?B=P.3*EPS5.MR4-R3) 

UB = U3*EPSij° (U4-U3) 

V8 = V3*EPS5<' (V4-V3) 

..SB = S3>EPS5» (S4-S3) : 

UC = 0 (NLC.H) ®COS Th£ (N rC)-V (NLC.M) oStNThe (NcC) 
VC = >J (MLC.M) *SIMhE <NrC)'*V (NLC»M) •CoSThe (NCC) 
-P_tN N A « M i. =P„( Nl C ♦ M)-«E p c f>«( PR-P-VNLC.M h ) 


— NOS01560 
NOS01570 
NoSO 1580 
-NOS01590 
N0SOI8OO 
NoSO 1610 
-NOS01620 
NOSO 1630 
NOSO 1640 
-NOSOV650 


S (nNa iMIsS t)NLC ♦ t H>SP-c6o-(-S8 ~ S-UVLC-+M4-) NoSO \ 660 

‘ UAsUc+EPSP" (UR-UC) NOS01670 

VAsVc+EPSt)" < VR-VC) NoSO 1680 

U-(nIN4 ,m) =UA«COSThE ( W^A>-^VA«SImTHF U‘AM- NOSO 1690 

45n V (NNA »M) = -L A»StNTHE (mAa) *VA*COSThE (NAai NoSO 1700 

49o CONTINUE NOSO 1710 

00-498 -M=l «MCfl 

00 4g 2 Na3tNCfl 


NNsNfiEG ( 8 ) *N 

ACH *SQRT-l tu_( fsns, — ,M > «*2-.V4lslN— » m-> a «3|V4GAMPa*E'XP-<P+Nn-,M)-/G A*S <NN-»M 

1 ) /GAPMA) ) ) 

IE (ACH.LT. 1 .01) GO T* 492 

Nj.aN_.V-- _ 1 . 

NL1=NN-1 „ . 

ACHlrSQRK ( 0 (NLl,M)«a 2 .v(NLl»M)«* 2 )/<GAHM 4 *ExP(P(NLl*M)/GA.S(NLl.M 

44 . / gamma.).). 1 

IE (AfHl .lT. 1 .01 ) GC rO 494 
492 CONTINUE 

60-4 G 499-- — 

494 N2s‘>2 

NL2=NN-2 - 

UC N- s U < N NT P. > « C 0 S T H£ ( u) _-V.( NN i M ) • S I N tHE4-N > 

VCNl sll ( NM « M ) *S INTHE (.,)*V(NN»M)*COStHE(N) 

UCn2 = U (NL2 »M) *cGsThE fN2) -V ( NL2 f M ) I NT hE (n2) 

___VCN2sU.(NL2«M)-*SlNThE.,N2) ♦V-(NL2TMvAr-OSTHE-(N2) 

P (NL 1 * M ) = < P ( NN « R) + P ( »'L2 » 8 ) ) /2 . 

S (NLl »M) s <S <NN»M> *S ( mL2 »M ) ) /2. 

_XlCNl = (UCN*LCN2) /?. 

VCNls (VCN+VCN2) /2. 


U (Mi.»M) = UCN1*C0STHE (N1 ) »VCN1*SInTmE (Nl ) 

_y-(N ; Ll ,M) =rL.CN.lisiNTHc4tvU4-4-VCNL*X^SrHE4NU 

49fl continue 
499 CON [ I NUE 

5.00— CON f I NUE 

C INTERPOLATION AT INTcRIOR POINTS 

00 'QC L = 1 * 2 

NST Aa2 

NEINsNMC (L) -1 
00 /OO NsNSTa.nEIN 

_____ NRs (L-l ) *N?'A-X.N 

IP (L.EU.2) 00 TO 520 
NNsNREG ( 7 ) *n 

MFlNsMf'C (NR* .1) Vl 


MSTAr l 

GO TO 530 
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-NOS01720 
NOSOI 730 
NOSOl 740 
-NOSO 1 750 
NOSO 1760 
NoSO 1770 
-NOSO 1780 
NoSO 1 790 
NOSOlaOO 
-NoSOI 810 
NOSO 1820 
NOSO 1830 



S2JL_PS.TA=MmC ( NflJ2l«-l 

WFlNsMCJ 
NNsNrEG ( 3) *N 
S3fl qQ . J n 0 - -M = M SJ 4UM E-t ti 

JJsOp (NR t M) 

LL1=L1 < JJ) 

LL3=L3 ( J JJ 

1 1 1 = 1 1 (JJ) 

I I2 C I2 < JJ) 

£6-S3_cEP 3( J Ji 


-N0SOV84O 
NQS01850 
NOSOl 860 
-NqS0 1870 
NOSO 1880 
NOS01R90 
-NQS01900 
NOSO 1910 
NOSO 1920 
4*OS01930 


.EPS 4 sEP_4 CJ JJ 

LJl=LLl-NHF(j(fl) 

LJ3=LL3-N«E(j (R) 

UlC=LI.(LLl-»IiU*C0SjHc--(LJ14^aLU^rr4-) Jt SIN-THF-<L-JJ-)- 

U2C = U<LLl»Il2)*COSTHr(LJl)-V(LLl.Ir2)«SINTHE(LJl) 
U3C=U (LL3* Hi )«COSTHr (LJ3)-V (LL3 , 1 r 1 ) *8 INTHF (L J3 > 

UA.C-=UtLL3!.ll2)*CqSlHe-tLJ3)J5V:.(LL3*tT2^*SINTHE-(LJ34-- 

V1C S U (LL1. till )«SlNTHr(L.Jl) *V <LL 1 , I T 1 ) *COSThE (LJ1> 
V2C S U <LL1 t U2) oSINTHc <LJ1 ) *V (LL 1 . I t2 ) *COS THE (L J 1 » 

V3c=lULL3»:m l^SlNTHr lL-J3r>-V (LL3^-I t 1 ) *COSTHE <1 J3 ) - 

V4C=U ( LL3 » 1 1 2 ) *S I NThr (L J3 ) ♦ V <LL3 . 1 T ? ) *COSTHE ( L J3 ) 
P8=P (LL 1 »II^)*EPS3«<r.(LLl«II2)-P(Li 1*112)) 

.U0S.O2C* LP S3?C.U4C-U2Cv - : 

VB = V;>c*EPS3* (V4C-V2C', 

S8 = ?(LL1 .ir2)* E Ps3*(c(LLl.n2)-S(Lj 1 * 1 12) ) 

__P_C=?4LU . I U-) *EPS3* < r UU.-3* J-I U-P-<4t l -»4 -1 1 -M 

UC = UlC*EPS3* (U3C-U1C', 

V.CaVlC*£PS3 # (V3C-VlC» 

SC»S (LL-1 « 1 I M * E P S 3 *-(-< + ti. 3»-M4-W5-f 

•P (Ni'J*M)s'P«*tps4* (PC-oB) 

S<NN«M)*Stt*tPS4*(SC*«0) 

U.(M*i t M)=uB*fc.ps4« (UC—tiB) — - — : 

V (NN»M)=VB*tPS4*(VC-i/B) 

700 COnUnUE 

R E T Qp M 

END 


-NOSO l 940 
NOSO l 950 
NoSO l 960 
-NOS01970 
NOSO 1980 
NqSO 1 990 
-NOS02000 
NOS02010 
NOS02020 
-NoS02o30 
NOS02040 

NOS02050 
-NOSO 2 060 
^OS 02070 
NOS02080 
-NOS 02090 
N0S02 100 
NOS021 10 
-NOS02120 
NOS02130 
NOS02140 
-NOS02 150 
N0S02160 
NOS02170 
-NOS02180 
NOS02 190 




SUBROUTINE ...KLUBO.(X» y'.yprj 

C0MMCN/HLK5/XP(200) ♦ CP ( 20 0 > * KE » L I Tp 
DIMENSION ft <200) * B ( 2 * 0 ) 

DIMENSION AW (20 ) 

OAT* LAMP/-/ 

OATA LI PE/ / 

IQ O-F OR M ft T-( 2 0 A A 1 : 

101 FORMaT ( iiT 1 *4) 

102 FORMAT ( 1X»3H,4E16. 4', 

JJQl3._F0RMaT..( 10 1= ) : 

IF (LAMP. EC. C) L I Tg =j IP£ 

IF(lAMP.EC.I.aNO.LItE.NE.I) GO TO 50 

IF I Li T E . EC . 1) -60—10- 2 5 

REaO(S.IOO) AW 
WRITE <6, IOC > AW 

XP-ll 1 = 1. 

■ YP ( 1 ) * « 5 
REAU ( 5 . 1 o3 ) KK 

W R I T F ( 6 *1 S 2 J — K K 

KK1=KK*1 
00 l(i K = 2*KK1 

— , RE AU-( 5 .1 0 1 ) XoP-»-YPP 

IF (EcF (5) ) 20,2 
2 CONTINUE 

xP-( K-i=xPP/2-. *-i , — — : . — 

YP(K)=YPP/2. 

WRI^E (6,101) XPP.YPP 

-Jul-xonlinue r : 

2o k£n u =kk 1 
GO TO 29 

_2S-KENQ=KE_ , 

L I TE = 2 

29 IF (KEN0.LT.3) go TC iO 

KENUlsKENOfll - 

KEM)p = KENU*l 

A ( 2 ) a ( ( YP (3) - YP (2) ) / 'f XP (3) - XP (2) ) - ( YP ( 1) - YP (2) ) 

1 _>l< * R < 1-1 - x P U )-)-) /.< X F < <* 1- x P ( M ) 

B(2>a(YP(l)-YP(2)-A (?>*(XP(1)-XP<2) > **2) / { XP ( 1 1 - XP ( 2 ) ) 

IF (KENQl ,LT *3) GO TO 35 

00.3(i_K = 3, KEmo 1 

0(K)r2.*A(K-l )«(XP(K', -XP(K-l) ) ♦B(K.l) 

30 A (K> a < ( YP(K*1 )-YP(K) \ -0(K)* ( XP ( K ♦ 1 ) -XP (K ) ) ) /<XP(K*1 ) 

1-XP-lK )-) *«2 

35 XP (KENOP) =XP (KENO) *1 f 

8 (K^ND) =2.«A (KENDl ) * fXP (KENO) -XP (K fNOI ) ) ♦R jKENDl > 

A(kEnO) =-B (^END) /-(2.*-(XP.<kEN0P1-xP-(KEN0-)-)-)- - 

XEnOsXP («,ENUP) 

YENlja a (K£NC ) * (XEND-Xo (kFnd) ) **2 + B (kEND) * (XEND-XP(KEND) ) *YP (KEftD) 

_WP 1 1 E-( 6 il 02J_(K iLAMR.ti.I.TE , xE-(K)-t.YP.( w_M A (K 1 , B-(K > ,K=2 ,KEND 

40 LAm*=1 

50 IF (^ENO.LT.3) GO TO oO 

I F.( X . G E • X E N 0 ) _g o -T C—5 5 

DO Of) K = 3»KtN0P 
IF (A.LT.XPJR) ) GO TO 7o 

JbQ_CONT INUE-— TOT : 


pluotoio 

PLU00020 

PLU00030 

PLU00040 

PLU00050 

PLU0O060 

-PLUOOOTO 

PLUO 1080 

PLU00090 

PlUOIIOO 

PlUOUIO 

PLU00120 

^LUOn 1 30 

PLU0O140 

PLU0C150 

-P LUO-1 160 - 

PLU0H70 

PLU00180 

PLUOi 1 90 

PLU00200 

PLU01210 

PLU00220 - 


— PLU04230- 
PLU0O240 
PLU00250 
--PLU0O260 - 
PLUO 0270 
PLU00280 
— PLU0O290- 
PLU00300 
PLU0O310 
— PkU0o320 - 
PLU00330 
PLU00340 
-PLU00350 - 
PLUO 1360 
PLU0O370 
— PLU0O380 
PLU0O390 
PLU01400 
— PLU00410 - 
PLU00420 
PLU00430 
— PtU00440— 
PLUO 0450 
PLU00460 
— PLU0O470 
PLU004B0 
PLU00490 
- PLU00500 
PLU0O510 
PLU00520 
-PLU00530- 



^O-LL**- 1 PlUO'1540- 

Y»A(LL)*(A-Xp( L L))««9*8(LL>’*(X-XP(j L))^YP(LL) PLU00550 

YPR*2.*A(LL>*(X-xP«Li ) >*8<LL> PLU01560 

— RE TURN — PLU00570 

75 Y=r^\D PLU00580 

YPR = 0. PlUO 0590 

RETURN PLU0O600 

80 Y=YP(1J PLU00610 

YPR s 0 • P|_U0 1620 

— — PE T URN p L U 0 n 6 3 o 

EN 0 P(_U0 1640 


SUafiCU.T-IfclE_OU.TP. T)ut0 

COM^CN/ULKT/R <;o*2) ti'i(3i.?l*W(3o,2i tP<3n»?) »S(3o«2) «PN(3o*2> OuToi 

1»Un(30»2) hN(3o«2) «Sm(3o» 2 ) » T (30,2) *PO(30«2) *UO(30«2) ♦ WO <3o * 2) CUTO^ 

2 ♦ SC t 3 il* 2 A., k A .gAVMa *GA .Rb»6C>6rv»«E *GP,Ga«Pl»TT0T GuT-010 

3»02*K.DX(2)*D0x<2> * NrM ♦ C ( P ) *0 ( 2 ) , C»' ( 2 ) , BN f 2 ) , I « I A tfcjA ( 2 ) • BZ { 2 ) OuTOiO 1 . 

4.C2‘2)«X(3 *2) «S0E*B!rTAP»?*NASAVF*7NtL1 »NSO*NS»KCUnT*OPMAXO OUTOI 16 

COMMCN2ULK2/PRATlO,PrvIST-«PRAO,KPl,-UuE» JJ QUT0'10 70 

1001 FORMAT (///AX, i^HpLUMr STRUCTURE /4 v ,5HXSTFP , 15 ♦ 5X t 2HX= OuTOlOSO 

l»F8.3«5X,3hUX=,El2.4 / /9X«l HY , 9X f 1 Ho ,9X * 1 h V .9X . 1 RU f 8X * 3HTAU 0 uTC v 0 90 

■ — 2 *8 x • l HM» 9X «.i HS *9 x 1 1 Ht } T>t>T0il00 

1002 FORMAT <4X»'jriR EG i 0 N «?4) OUTOll lO 

1003 FORMAT (13* HF1 O.S) 0UT01120 

_ WRITE. (biLOTU-K^Z, 02 0UT01130 

00 A 0 I = 1 • I A 0UT01140 

WR I TE < 6 * 1 0 12) t 0UTOO150 

hlCsMA T4-MO- — OtfTO >160 

IF ( I .EO. 1 ) fc.C=NC*l • OuT 0 H70 

00 lo U=l*fcC OuTOilSO 

IF-f-i «N . EQ» 1-lGO-^O— TO OyT On 1 90 

TAU=U <N. I ) /* (N, I ) OuTO 1200 

AMaCh = S'JRT t <U (n » I ) ««?*W (N , I ) **2) /GaMMA/T (Mt I ) ) OuT 012 1 0 

P R E S = F-x P < P < N , . i y ) q U T 0 1 2 2 0 

WRITE <6« 1 0 • J) N,R <N. I', tPRESiU (N» I ) iu <N? I> » TAU* AMACH , S (N « I ) « T (N« I ) CUT 0^230 

In CONTINUE OuTOn240 

— RETURN . ^ : OuTOi?50 

E r, 0 OuT 01260 
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vt uu r i ii iw,w*r jvc uriai 


U3/V1/ f J 


SUBROUTINE p LHME PLU0O010 

C PLU*E OEC* ' PLU0O020 

C0mMCN/«L«7 / R (30.2) *u<30*2)*W(30.2mP(30*P)*S(30«2)*PN(3o»?) PlUOoo30 

VtUM30«2M*Nr3 0 »B) .S». <30*2 >*T< 30.2 >-»P0<30.2>.U0< 30,2 > .tf0(3o»2> — - PlUOtq^O 

2 ♦ SO < 3 0 * 2 ) *^A,KA,GA^M^,Ga*r,ti.GC.Gn»rtE*GF*6O,PI*TT0T PLUOnoSO 

3,DZ*K.OX(2>*DOX(2> .NrM.C ( 2 ) *8 (2) .Cm (2) ,BN J2) . I *IA*nA <2) *PZ (2) PLU0T060 

4.CZ l 2> * * <3 *2) .SOP*BrTAF.Z*NASAVF*-7N*L-V*NSO*NSrKOUMT*DPMAXO PlUOOo70 

C0mMCN/oL*F''PpATI0,Pp.I<;T,PRA0.KPlUwE, Jj PtUOnoBO 

CO(-8cN/HL^5 / ZP(2oO) *rP(200) *«P*LITp PL<J01090 

DIMENSION- PM 20 Q^*P*f2oO>*OK-< 200 >-*wK-t2 ; TO0 PLU0O1 10 

NAm§l IST/JE T /nA ,mA ,K a .KOUT.PRATIO. nlST .GA mmA .STAB .aCH.TTCT.KMAP PlUOO 1 00 
0AT“ LAMP/ / PLU0nl2o 

— M>01- FORMAT (///AA-, V6HPL~l‘Mp -INPyT— OA^A-A/4-X-.?HPr>ATIO=*F6,3.5x*5HMACH = PlUOt 1 30 

1 *F6,3 «SX,5h fTOTs.FP . ->//4X ,6HNA ( 1 ) =, I2.5X.6HNA (2) = * 12 * 5X . 2 HmA = PlUO'1 140 

2» 12 *c X ,3HKA=. 14.5 X.5uK0'JT=. 14 ,5X.SuKPAP=, t 4/4 X ,5H0i ST =. F t .3 PLUOo 150 

3.5X .6HGAMM A s.p*,, 3*5 Xj-SHS7 ABs-, P6.-3) PLU0M60 

1 002 FORMAT (///4X, ! 4HpLLMr BOUNDARY // 1 i X . 1 HX . QX . I H Y ) Pl u 0H70 

1003 FORMAT ( 1 X» IA,2F10.5) PLUOolBO 

IF-fLAMa.NE.-U4— 60— IC— f 0 Pt_U0<i 1 90 

NA<l)a2i) PLU00200 

NA < 2 ) a 1 5 PLU0T210 

MAbIq P L-UO ^220 

KAsloO ' ' ' PlUO 1230 

0 I S 1 = S . PlUO'IBAO 

.GAMMA =-J-*4 PLUOO 250 

STA8= 1 . PlDO o 260 

ACH=2. PLU00270 

KMaP-=«1 P LU 0^280 

PEaD(5.JET) PLUCT290 

wRI]f <6. 10 .1 ) PRaTIO - , ACH.TT OT.N A, Ma.KA, KOI JT.KMAP.DIST, GAMMA, STAR PlDO' 1300 
L-AM-P-s 1 Pf UO 0 3 TO 

NA1=NA(I) 

NA2=M (2) 

Lfl_C 0 Nil N U E — Pf-U 0-T320 


NA{1)=MA1 

NA(2)=NA2 

-Lite = 1- 

LL = 1 
NS = 3 

-NS0 = 3 - 

OPmAxOs.i, 

kou*t=c 

_2P.flVsP.DlSI 

CP ( 1 ) sPRAU 

NAS?vE=nA (2 ) 

-GBsl ,/ (GAMMA. v> ) _ 

GA=GAMMA*Gr! 

G0=;5/GB 

_GE a 1,*GD 

GC=&E/GO 

GFsScPT (GAMMft) 

-GGsScRT < G C ) 

PI=4,«AT AN ( l . ) 

ZN=Pn 1ST 

CAL-L -BOUND fZN ♦ POUM V 


— PL-U00330 

pluo'vsao 
PlUO^SO 
—PL-U0I1360 
PLU0-1370 
PLU003B0 
— PLU00390 
PLU00400 
PLU0T410 
-PLU0O420 
PLU01430 
PL U0T440 
— PLU0T450 
PLU0'i460 
PLD0T470 
-PtUO-iABO 
PLU00490 
PLU00500 
-PL-U0^510 



- - <VV • I'l ur 1*1 


ya/v»/ f J 


PAaP CiiM- 

TTOTsl .♦GU«ACh«*2 
POsTtOT**GA 

SMF“(TTOT«E^Pt-PA/GA',^U^/GD 1 

DZa«o5 

OUM* (SMF-ACH**2) /2oO* 

00 -2n_K*l t2U& 1 1 

00M=K*0UM 
DEP a S0RT (OCM) 

SMK = ACH.**2.+JJQm : 

CU=ACH**2-i. 

C22=SMK-l. 

ENUN=GG»A-T-AN-(-<;gRT4C1 ? /GC-> 7-A7 AN4 sQo44CAV44 

ENU^aGC^ATAN ( sort (C 2;/GC) ) -AT AN ( sQoT <C22) J-ENUN 
PSIK = Enuk-a Tan ( 1 ,/SGoT (C22) ) 

OUf.=4-«-+.GO*£MK 

RK (* J =PWAO*UZ*TAW (P-StK) 

PK ( K j sGAe-ALGG (TTOT/Dr^) 

— OKb^GR-T-< GAV-MAoGE^SMK 

WK (N) =OK*CCS(enijk) 

UK (Kj sOK*SIN(ENUK) 

2.a_C0vj-iNUE 

RFsPfiAD»OZ« T an (ENUk j 
ORsRf/NA ( 1 

NCm=nC-1 

00x(D=na(1) 

DA(J J = l*/UCA t u 

00 8 0 N=2»NC 
R (N» 1 ) = (N-2>*n& 

JLtN»44.a-(N--21 *ox < 1 ) 

I F < H ( n . 1 ) .GT.RK(l))Go TO 30 
K = 1 

EPS = 0. , 

UK < * ) * •* . 

WK (Kj =ACH«GF 

PK( K )_aiu 

GO ]C 70 

3o IF (R (N* 1 ) .LT.rk ( 200 > \ Go TO AO 

Ks.J ’9- 1 : 

EPS=1. 

GO Tc 70 

AO — 0 0__S o- K_= U-LV-9 : 

IF I (r(N»1)-Hk«k) )*CR(N, 1)-RK(K-»1) ) ,LE.n.)G0 TO 60 
5o CONTINUE 

_6j 1_EPS= {P_.Ifct*-1 1 ”RK (KJ 1 /.(oK ( K ♦! 1-RK1 K4-J 

7 0 U (N* 1 ) =UK (K) ♦EES* (UKfK*l ) -UK (K) ) 

W (N» l> =wk <K J *EPS* O.K jK«I J -WK (K) ) 

PIN »JJ_8 PKIkL+EPS AtPJlt-tU 1 ) -PtCtfUJ 

S ( N * 1 ) = 0 « 

X (Ml ) =UX ( l ) * (N-2) 

. 6 0 — T (N »1 J.aEXP (P-tN* 14/6A*S /-GAKMA-i— 

Ull'l )=-U(3*l) 

W ( 1 • 1 ) s* <3» 1) 

S-U » U = S43 . 1.) 


-PtU0oS20 
PLU0OS30 
P|_U00540 
-PLU0O550 
PlUO' 1560 
PLGO'1570 
-PLUO'tSBO 
PlGOOSRO 
PlUO^GOO 
-P 4 .U 006 IO 
PLU00620 
PlGO 1630 
~PtrUoo640 
PlUO 0650 
PL^O 0660 
-PtG0ft67(T 
PL^O 0680 
PlGO 0690 
PLU0O700 
PLU00710 
PLU00720 
_4>tU0o730 
PLU0O7A0 
PLU007S0 
^^0076(7 
PLG00770 
PL^O 0780 
-PLU0A790 
P|_U0"800 

plug > a i o 

-PtU0o820 
PLU00830 
P|_U0o840 
-PirU0ofl5O 
PLG01860 
P(_UO 0870 
43 14)00880 
P|_GO 0890 
P|_ u 0o900 
Pt4J0O9l0 
PLUO 0920 
PLU00930 
^t:U0o940 
PLU0OR50 
PLU0O960 
PEU0097 o 
PLG0O980 
PLGO 0990 
-PLUOlOOO 
PLUOlOlO 
PLU01020 
PtG0l030 
PUUOIOAO 
PlUOI 050 
PtU01060 



uu r hi »y,u-rJWC U*" I 3 l U^/US/73 


P I-)., i ) =p ( 3 ♦ 1 > 

T(l»l)=T(3»l) 

BETAp=SURT (S mp-1 .) 

SQp= -smf«Gamma«T (NC » v) 

U ( 2 ♦ l > = ■ C . 

K = 0 

Z=D^*PDIST 

I As 1 
8 ( I ) sO . 

U '(-*-»■}••) =-U 1 3 »}} 

X ( 1 * 1 ) =-0* ( 1 > 

8Z(i)si. . 

€-M-Ib94N0»4-1 — 

CZ(l)=<C(i)-P9AO|/OZ 
CALL outp 

L 1 = U : 

100 K=K*1 

ll=ll ♦ V 

OZs'i-, 

00 lit) 1 = 1. 1 A 
NC=NA(I)»1 

IF~<£ ,EQ • 1 1\C = mc*4— 

oo lir ns 1 »nc. 

DEN = w (N« I ) «*2-GAMMa#T (N, I ) 

DUM = SORT (.IL-IN, I 1 «*2 *. rf lM» I L«A24/<»AW*»A/4-m» T-WI-t-I 

0 1 p s ( - hZ (I)”X(K« l)*<rZ(T)«BZ(I) ) ) 

OEw = flRS (01 '* <11 (N, I ) <N, I ) *G ANM A*T ( N » I ) «0l IM ) /OEM 

DAV = ABS-I0i^-UMN»1) M--sGA«M4«-7-<N* I >-*DuM)-/OEM 

IP (DAM.GT.tt-M) 0EM = 0 a M 
OZl=STAb*ux ( I) « <C < l 1 _B ( in /OEM 

IF (021 .GT GO— TO-OvO 

07=021 

lln CONTINUE • 

1 F- ( * 2 . 6 L, 1 „ E =A LGO-IO-12-g 

CAl l OUTP 
STOP 

J-2/i--ZN=^*0Z : 

CALL SUPEH 
ZP(Ll)=ZN 

cp-(Ll) =c (I a 1 : 

130 continue 

IF ( ’Ik/KOijT) "KOuT.EC.,?) call outp 

IF 1 1 .LJ ..OllF. and .K.Lp.KA J-GO-XO-VOT; 

KPsLL 

WRlTe(6tlO,; <?) 

1 t>-L Jl (KK*2P-(i?K)-*CP-tKK-L*K«sT-»KP-i : 

RETURN 

ENO 


-PLU01070 

PLO01O80 

PLO0109O 

ploouoo 
PlOOI 1 10 
PLO0112O 
-PlUOI 130 
PLOO 1140 
PLU01 150 
-PtOot 160 
PlOOI 170 
PlUOI 180 
-PLUO 1190 
PLUOI200 
PLU01210 
-PtL»01220 
PLU01230 
PLUOI 240 
-PtrUo 1 250 
PLU01260 
PlUO 1270 
PLUOI 280 
PLU01290 
PLUO 1 300 
-PLUOI 310 
PLUO 1 320 
PLUO 1 330 
PLU01340 
PLUOI 350 
PLUO 1360 
-PLUOI 3 70 
PLUO l 380 
PLUOI 390 
-PLU01400 
PLU01410 
PLUO 1420 
-PLU01430 
P.LUOl 440 
PLUO 1 450 
PLUO l 460 
PLUO 1490 
PLUO 1 500 
PLUO 1510 
PLUO 1520 
PL u O l 530 
PLUOI 540 
PLUO l 550 
PLUO 1560 


125 




_sub Routine super 


— — SuPOnolO 


COMf'eN/HL'W/'R (30*2) *M <30.2) *W< 30 . 2 ).P< 30 »?).S( 30 » 2 )*PN< 30 » 2 > 
l *UN <30.2) » (3o«2) .S.' (30.2) «T <30 .2i ,PO (30,2) ,U 0 (30,2) ,WO (3o.2) 

2 . SQJ.30-.-2l * V * iKA tGAMMn .GA,Gfi-, 6 C .G(l.A£»QE-» 6 fi,P-I » TT-O-T 

3. DZ.K »Ox (2) *00x ( 2 > »NrM,0 <2) » B ( 2 ) , Cm < 2 ) • BN ( 2 ) , I , j a «nA ( 2 ) * 8 Z ( 2 ) 

4, CZ ( 2)»*<3 *2) »S0F«8cTaF.Z»NASAVE »7N*Ll » NSO. NS . kOUnT . 0PM A XO 

COWMCN/ bLKE/P-RA-T 1 0 . Pr, I PR AD^KPtUuE » JU : 

DIMENSION oh (3o) ,02P ',4o) ,C0 <2) . ToY<2) .ERR (2) 

100) FORMAT < IX. 3UHWARNING’. I TEPA T ION FAjLS AT K= . IS ♦ 2X ♦ 2 HZ= .F 1 0 .5 ) 

_1 002 -f 0RM4T-(2 qI£) i 

PB0Un=P(NC.1A) 

IF ( JJ , £D. D ) GO TO 60 

c LF(jLl. £ Q >1 . 1 -GQ-XQ-fcO— 1 

0 PwAX=.w 
NC=Nfl ( 1 ) *2 

00 Lf) Ns.4.* \C : 

PGR=P ( N , 1 ) -P(N- 1 ,1 ) 

IF (H g r.LT.CPmaX) GC TO 10 

DP M A-X--S P G S : 

NS = N 

If CONTINUE 

1 E4AS^ME. NS (D_K 0UNT = * 

IF(Nc.NE.NSO) go TC *,0 
IF (0PMAX.LT .DPmAXO) «0 TO 60 

kountskoijni-* i---- _ : — . — 

IF (fCUNT.LT ,3)GO TC *0 
20 CON I INUE 

CAL L - OU-IP- : 

IAsf 

NCSjsNA < 1 ) *2 

N stvS->3 '■ 

NCsNA ( 1) .2 
OCX ( 1 ) =NA ( 1 ) 

_OX-(-Li-= L. /UOX 44 -) — 

C < 2 > rC < T) 
cz (2 > =CZ ( 1 1 

__ — — C(L)=P<NC»D 

CZ ( l ) = < U (NC * 1 ) «W (NC *7 ) -GAMMA*T (NC» i ) • SORT ((U(NC.1>**2+W(NC.1)««2 
I ) /GAPMA/T <-'C» 1 ) - l , ) ) / (w(NC» 1 > «* 2 -GaMMA*T (NO* 1 ) ) 

8 Z-< 2) =CZ <11— : 

B ( 2 ) =C < 1 ) 

NC»'^A(2)*1 ‘ 

OOx -<2 ) »na (2 ) — : 

OX (C) = ) ,/U.QX ( 2 ) 

D 0 4 0 N = 1 * •“» C 

XL.N..-2 ) -OX IEL<UN-LL i 

R(N.j)=«(2)»(C(2)-B(^))*X(N.2) 

NSm s NS-1 

... I P P -0 = (1 : : 

00 35 NNsN^M^cSA 
DUM S R <N»2) -K <NN» 1 > 

I F_( iPPO.EU ,.U_ G 0 — IC — 3 0 

IF (^LM.GT. . 'J) GO TC -,0 

EPs (R (N,2) -R (NN-1 . 1 ) ) / (R (NN. 1 ) -R (Nm-1 ,1 ) ) 

P4n* 2)-=P<N.\-1 ,^-).*EP«>P(NXf^l-)-P-<LlM-V,lH — — 
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Su p 0 ip 2 o 
SuPO'10 30 
-SuPOOOAO 
SuPOf 050 
Su p 0l060 
-SuP0no70 
SuPOOOBO 
Su p 0l090 
-suPonoo 
Su p 0 H 10 
SuPOO 1 20 
-SUP00130 
SuPOo 140 
Su p 0n 1 50 
-i SuP-0 •» 160 
SuPOf 170 
SljPOilBO 
SuPOf 190 
Su p 0 l 200 
SUP00210 
-SUP0O220 
SuPOf 230 
SUP0O24O 
-SUPOO 250 
SuP0f?60 
SUP00270 
-SUP0O280- 
Su p 0O29o 

SuPOf 300 
-SUP0O31O 
Su p 0l320 
Su p 0l330 
-SUPO0340 
Su p 0i350 

SuPOo 360 
-Sll p 0l370 
SUP00380 
SuPO 1390 
-SUP0O400 
SU p Oi4lO 
Su p 0i420 
Su p 0 ; '430 
Su p O 044 0 
Su p O 0450 
-SU p 0-1460 
Su p O 1470 
SUP00480 
SU p 0O49O 
Su p 01500 
SUP0O510 
~SU P 01520- 
Su p 0i530 
Su p 0i540 
-Su p 01550 



U ( N ♦ 2 > =U < NN- -1 1 V 1 -*EP 14-U < NN-v* 111 

tat (N*2 1 = * <NN-1 , 1 ) *£P* ,W (Mf». 1 ) -H (NN-1 « 1 ) ) 

S (N»2> =S (N\-l ,i ) *£P* /S fNN.l ) -S (NM-1 ♦ 1 ) ) 

TIN # 2 1 =ExP-(P4N,2>/GA*SfN.2E/6AMMA> 

IPPOsl 
3 0 CONTINUE 

-40—CON 1 I NWE— 

NCs^A (1 ) *2 
00 S() N=2»NC 

X-{ N t H s' < N-2 i 

50 CON T I Nut 

P(l*2)=P(N?M,i ) 

U4F»2 y a O' ( MSM.-H 

V»(l»2)aW(NSM,1 ) 

S(l»2)=S(NSM.1 ) 

T44 -» 2 ) -* T < N SM»n 

NCsNA ( 2 ) ♦ 1 
P(nC»2)=P('.CSA.1) 

U4nC .2 > =U t N. Csa *-1-) 

W(NC,2)=W(\CSA.l) 

S(nC,2) =S(NCSA,1) 

T- INC. 2-)-=T- in c s a .1) 

Llrl 

CALL, outp 

-6i_ CONTINUE — T 

NSO=NS 

OPm ^ xOauPM A A 

— : LO C P-= k- 

7n 00 2\ 0 1 = 1 * IA 

nc= n a(I)*i 

: IF-4-1. EQ .4 LivC=MC*4 

00 ton N = i.A/c 
IF(lN*I.EQ.1 ) GO TO 2*0 

NM.J =N«LOOP. : : 

IE I I s ! • EO • 1 ) “> M 1 =N 
if(N,eq.nC)Nmi=n-i 

NP 1 =NMl M 

PX = <p INPH 1 1 -p (NmI *Ii ) *nox ( I 1 
UXs «U (NP1 • I > -II (NMl . I ) ) *00 X ( I ) 

WX.=_Lw.( .1 *ncx_LIJ 

SX= (S (NP1 • I > -s (NMl . I'l ) *O0X ( I ) 

GX=*X/w(N.I» 

A Aa4_. /.IM. ( N.i. 1 ).« * 2-G AMn A* T ( N • 1 14- 

. IF (IA.E«.2.«NI).I.EG.7.ANn.N.EO»NC) GO TO 140 
I F ( I .EQ.2 .and.n.eO.Na (?) ) SX= (S(N, 2i-S (N-l ,2) )*D0X (2) 

IF_(J-.NE-«1 • CM *tJ *NE-»2 ) r.Q — T-O-B-O 

PX.= 0. 

sx=o, 

XX: — 

QX = C. 

8 A TAu = L'(^*Il/ w (N.T) 

A8 = AA«'a.(Nt.[.) **?*-T-AU 

AC = Gamma*ma*w (k , j j 
A 0 = AA*.TAIj«T (N, I ) 

xr= !*/_( ci lira cU-> : 


— Su p 00560 
SuPOoSTO 
suPooseo 
— Su p 0 o 590 
SUP0O600 
SUP 0 O 610 
— suPo^ftao- 
Su p 0O630 
SuP 0"640 
— Su p O ^650 
Su p 07660 
S(J P 00670 
— SUP00680 
Su p 0l690 
Sll p 0'>700 
— SUPO'1710 
SUP0O72O 
SuP 0 n 730 
— 5UP00740 
SUP0O750 
Su p 07760 
— SUP0-A77O 
Su p 0o780 
SU p 07790 

— Su p 0o800 
Su p 07810 
SU p O 0820 
— SU p 0oR30 
S(J p O 0 84 0 
Su p 0O850 
— Su p 00860 
Su p Oo07O 
SU P 00880 
— Su p 0 >890 

su p ooooo 

Su p 0G9lO 
— $U P 00920 
SU p 0O930 
SuPOnPAO 
— SU p 00950 
SU p 0'i96O 
SU p 0A970 
— SUPO^NSO 
SU p 07990 
SlJ p OlOOO 
— su p oioio 
Su p 01 020 
SU P 01030 
-Su p 0l040 
Su p O l 050 
SU P 0 1 060 
— SU p 01o70 
Su p O 1080 
Su p 0l090 
— SuPOllOO 
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XZa An* IBZ ( Ii* (X-CW»l)—i*)*C2 U4*XtN.l44 

AE=X2>Ab*X= 

AFsX7*TAU*xH 

AG=GAMMA#Aa?X*N * J4 /W f N. t4 

IF(I.E0.2.AN0.N,eQ.1',G0 TO 100 
IF(I.EO.IA.Ano.M.EG.mC) GO TO 15o 

DUM=-AC«UXAAR 

DOF = AG<>UX*XK 

PZ3-aE«PX*GAmma*aB«Q;*XP*O0M*JJ*oU» 

QZ^AE«iJX>AO*Px»X.R*OoM*Xj«Li:OM 

IF (Jj,EO.0)«0 TO 90 
IF(Jj*I.E«.l.AKD.N.En.2)Gn TO 90 

PZ =H-2 »IA B * G A M M a /p( N fi-*nu 

QZsQZ* < AD*GAMma/R (N * r > -00M) *JJ 
9r SZs“AF*SX 

U-Za- A F AU X-XLN,-I-T/W-( nVI > O-X-P «PX 

WZ=*(N*I)*'Z 
GO Tc 1 BO 

-WtMF (LCOP.EO . l )-6»A-TG-^<^e 

MC1=NA ( 1 ) *2 

BETA l =SQRT ( <U( l .2) •*?♦«< 1 ,2) «*2 ) /GaMBA/T ( 1 ,2)-l.) 

BETA? =Su R-T- ( -( U4 2-» 2 VM 2-«-24-*-«-24-/G A M MA-/-T4-p-r2-F» 1-.4 

A A 1 = A A 
TAUlsTAU 

TAlj-2=U<2 ♦ZJZ.W < ?., 24 — 

A A2= 1 .✓(W(Z*2)«*2-GA m MA«T(2»2) ) 

ALAMl a (U ( 1 .2) «U ( 1 *2) _Gamma«T < 1 .2) *oETAi ) *AA1 

ALAM2= (U (2.^) *w <2.2 » .GAMMA«T42-»2> *a£ TA2MAA2 

EPS= (CM < 1 ) -G ( i ) -ALAM7*D7) / (OX (2 )/Xd. (ALAH 2-ALAM1 ) *DZ) 
IF (tps.LT. ,-M E P S = . (le 

PST AfisP.U.i c ) .* FPS« (P (,•». 24-P U . 24-V 

TAU^sTAUl TAU2-TAU1 ) 

AL=;5 

T AU-3 3 U C3.« 24 /.W. (_3_« 2 ) — 

PP7Ap= 2.*P <<2»2) -P (3.o) ♦EPS* (P(3«2).P(2*2M 
. TT(jSr2.*TAL 2-tAU3 + EPc*(TAu3”TAU2) 

P_S I A As. ALAES J. A R *. (X. -A i .) * P P..TAR 

TAuSsAL«TALii4 M.-AL)«TTiJS 
OUPlaW (1.2) **2 /BeT A 1 /T (1,2) 

D UM2= W ( 2. 2 )./ B E-T A2/-T-L? . 24 

DUP=0UMl .EPS* (0 Um2-DmM1 ) 

OOPl=!J( 1 ,2)* (l)(l ,2) (1 ,2)«BETA1 ) « a A l/BE.T A 1 

OX)P 2 r U L2-. 2 ) * UJ < 2 . 24 *m-( 2-. 24 ^ ETA-2 ) -« A A 2VAE-T-A2 

DOM=UOMl .EPS* (UOM2-OriM 1 ) 

DOM=DOM*GAVMA/c (1 ) «t) 7 

OA K = PS1 A R ~ lUTM *tT- A U S *0 aJYI 

KIP=1 
MEs A 

IftY (1-ls.C ZL1 ) 

TRy<?)=CZ< 1)*.S9 
Uo sqr=sort ( l .♦try <mE>«42) 

EMI *4 »/SQR : : 

EN3 s -TRY (ME»*EN l 
TA1=-EN3 

IA3=F-NJ : : 


— SuPOU 10 
SUP01120 
SUPOU 30 
— SUP0U40 
SUP01 150 
SUP01160 
--SuPOU 70 
SUPOUBO 
SUP0U90 
— SUP01200 
SUP01210 
SUP01220 
— SuPO 1 230 
SuPO 1 240 
SUP01250 
—SuPO) 260 
SuPO 1270 
SuPO 1280 
— SuPO 1 290 
SuPO 1300 
SUP01310 
— SuPO 1 320 
SuPO 1 330 
SuPO l 340 
-SuPO 1350 
SuPO 1360 
SUPO 1 370 
-SuPOl 380 
SU p 0 1 390 
SUP0U00 
-SijPOUIO 
SuPO I 420 
SUP0U30 
— SU p 01440 
SUPOl 450 
SUPO 1460 
-SUPO 1470 
SUP01480 
SUPO I 490 
-SUPO 1500 
SUPO 1510 
SUPO I 520 
-SUP01530 
SUPO I 540 
SuPO 1 550 
-SUPO) 560 
SuP0)570 
SUPO 1580 
-SUPO l 590 
SUPO 1600 
SuPOl 610 
-SUPO 1620 
SuPOl 630 
SUPO 1640 
-SyPOl 650 
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VUnCsUN CNC-1-* 1 ) *T4 1 *WM (NG1-*1 )-«TA3 

UWNCaUN (NCI * 1 ) «En1***m (NCI ♦ 1 ) *EN3 
SMnCsU*NC»°^/^AMmA/T ,NC1 « l ) 

PN 2 NsPN(NC 1 *1 ) ♦AL 0 G< 7 GA“MA‘»SMNC--fl 0 i-/Ge> 

UWN=l' w ^C* (GU* runic* l /GE/S p NC 
UN { 1 ,2) =U*» f -*EMl *vWnCATA1 

WN< 1«2>=U4N*EN3*VWmG»TA4 

SN(l#2)=SN(Nci,l)*PN'5N-PN<NCl«l)*GAMMA*AinG<UWNC/ , -IWN) 

TAU^NsUN ( 1 /toN{ 1 .2, 

PM < 1 ,2>=DA-' tDllM»TAU2v t 

ERR <wE> =PN ( 1 .?) -PN2N 
IF(ME.E(j.2) Go TO 12* 

V£ -f : : : : 

GO TC HO 

12? IFMPS(ERM(ME, , , L T, j E-6) GO TO J 3* 

: — I.F.4 A g S ( t R R - 1 d ). - £ R p U4 y . GT-*-l-« E- 8-) — GO — T-0 — T22 

TRY«=.R9*lR V (2) 

GO fc 124 

l22_TRVtJsTRX<-lJ--EPB4l4-«iiRY42)-J*-Y-4-l-TlV-iePR<2l--€RR< 1 M — 

124 ERR < 1 ) =ERR (^) 

TRY ( 1 ) = TPY (2) 

TRY. I p.J.a.T RY.A ! : 

KIP = K IP* 1 

IE ( K I P . L T « 2 *j ) GO TC TlO 

4R I-T E ( 6 « 1-6 ’ 1 ) -K ♦ Z— 

l3o CON f INUE 
CZN=TRYA • 

S0-T.C-2--O — 

IV; CONTINUE 

IETLCOP.E-’.l ) GO To ?00 

-CN.U.J =C-CU-*.CZ-(-M-«D2 

. NC = *A.< 1 ) *2 
NCm=nC«1 . 

__AA1 = 1 ./_( 4_L\C_. l.).ft*2-/ ? AMMA#T-l.NC iVLi 

AA2 = 1 ./ ( ti (NGM. 1 ) **2-<".Am m A*T (NCM* V) > 

BET A i =S()hT (lu (\ C »t)**2*W(NC » 1 ) <**:>) /GAMMa/T (NC 

— BET A2sSQRT.IJ UTNCM*.L>.**2*W ( NCH«4 l «*p 1-4GAMM a /4 (NC^ ♦ 1 ) - 1-. > 

A|_AM 1 P= (U<\G «]')«W(Nr « 1 ) *GAmmA*T (MC t 1 ) «BETA1 ) *AAl 
ALAMlMs C U .< N <r «l)*W<Nr « 1) -G A MM A *T ( mC ♦ 1 ) «PE TA 1 ) *A A 1 

_^ALAM2P.s (U (NCm, l ).*W INrM, 1) *G Amma*T-<»'C^ » l ) «RE-TA2 ) * AA2 

ALAM2M=(U(NCM,l)«W(NrM,T|-GAMMA*T{MC? , »l) <, BETA2)flAA2 
EPSPs <CN ( 1 ) -R tvCM* 1 ) _A|_AM2P*0Z) / <Dy ( 1) /XR* ( A|_'AM 1 P-ALAM2P ) *qZ ) 
__ EPSM= £CN.Ul-r:R{NCM.iJ_J m A4_AM2K*02J_/-<0-v-M4/XR*-(-ALAMlW-ALAM2Mj.#0Z> 

DU^ A =RETAl»T (Nr,i) 

0Uv2 = peta2<» T in cm, 1 ) 

_0.1 .= U CM2 ♦ E P. S P * (nUMl-0uM2) __ 

El=” (0UM2 *EPSm* (0UVJ„0U W 2) >. 

00n 1 =w (NC* 1 ) **2 

— D0m2-W-{.ncm* 1 ) o»2 

02 = UrM?*EPSP* fplMl-PoMz) 

E2=6cb2*e p 2M« (00m1-0^M2) 

_ P 5 T A R P sR ( NC M « l ) *EPSP*.(P (NC-»1> -P-<nCu* U-) : 

PST Arm=P (NC m * 1 ) *EPSM# (P (NC» 1 ) -P (NCm, I ) ) 

TAul=IJ (NC* 1 > /to (NC* 1 ) 

—TAu2 = ll(NCM t 1 ) /w (NCv ) 
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SuPO 1 660 

SUP01670 
SUPO 1600 

SUP01690 

SUP01700 

siipouio 

SUP01720 

SuPO 1730 
SuPO 1 740 

SUP01750 

SUPO 1760 
SyPO 1770 

S(/Po 1 780 

SuPO 1790 
suPoifioa 

SuPOlAlO 

Su p O 1 820 
Su p 0 1 830 

S‘JP0184 0 

S'JPO 1 850 
SUP01860 

SuPO 1 870 

SuPO 1 880 
SUP01890 

SuPOlPOO 

SUP01910 
SuPO 1 920 

SUP01930 

SUPO 1940 
SUP01950 

SyPO 1 960- 

SuPO 1970 
SUP01980 

SuPOl 990 

SUP02000 

SUP02010 

SUP02020 

SUP02030 

SUPO20AO 

SUP02050 

SUP02060 
SUP02070 
— — SUP02080 
SUP02090 
suPoaioo 

Su p 02 1 1 0 

SIIP02120 

SUP02130 

— suPoa iao 

SUP02150 
Su p 0?l 60 

SUP02170 

Su p 0 2 180 
SUP02190 
SUPO220O 


I A U S P = T.AU 2 * t p S p * t T. A U t- T AU 2 £_ 


SUP02210 


TAuSp = TAU2*tp<5v* ( T«UT-TAU2) SU p 0222O 

DEmIs- (U (NC M ) *R AMV A « T (MC. 1 ) /R (NC»T)« (IJ(NC*1 ) ♦W(NC«1 > «BET A 1 ) • AA1 i Su p 02230 
_OEm<J = -.(U.CNC^.U *6AVMv*T4MCKt l4VJi^NrM^-)-* t<MNCM^l-MW(NCM»l >«BETA2^- SlJ p 0 224n 
1*AA4> Su p 0225o 

.0|m 1»-(U(NC»1) «GAMPA«T (NJC»l)/R(NC*i)*<U(Nr«l)-W(NC«l) *RETA1 ) *AA1) Su p 02260 
-O IV<2s- (U<NCM, |> *gAVMa «T (KC**±) /« <NcM«-l4 • MHNC* rl >-"W <NCM* 1 ) •BETA2) - Su p 02270 


1«AA2) 

AKPsOEM2*EPi>P»(nEMt-nEM2) 

._AKK=0 IM2*EP-SM« M2 ) 

AKP=AKP«OZ 

AKM=AKM*DZ 

—US =akp»di «pSta b eu d 2«t au s p 

E3=AkM+E l* c STARM*E2«fAU5B 
TAuNs (E3-E1/D1 «D3) / (c2-El/01*02) 

-PALtNC »-U-= IC.3 - Q2«TA UMV/ 04 

TEr'PsEXP (P\ <NC* 1 ) /64) 

Q2=2. *GA« (TTot.TEHP) 

_WM (NC ,4 ) =SCKT <G2/44-**T4MN«*244 

UN <NC» I > =WK <NC»1)*TA,,N 


Su p 02280 
Su p 02290 
Su p 02300 
SUP02310 
SUP02320 
SUP02330 
Su p 0234O 
SljPO 2 350 
-Su p 02360 
SUP02370 
SuPO 2380 
SIJP02390 
Su p 02400 


SN(Nc«l )=••"■ 

GO_Io_2 CO 

l5o conT tnue 

call ROUND (ZN.PB8) 

P-H-( N»I4=P dfi 1 : 

PZ80UN = (PN (N, t j -pBCUvi) /02 
TAu*=UX/W If, ♦ I ) -TAUOQv 

SKFiLt XXOX* £ AP.t-.Pf4 (M i JV-GAl -Ui/Go 

T <M» l) =EXP(HN(n* I)/Ga*S(N. I) /GAMMA* 

BETAFsS0RT{Smf-1.) 

__S0FS$MF-ftGAMMA*;-<N*-I4 : 

alam = ae*ag«w (n, i) * re;af*xr 

TAuZs- (T (N, 1) *pFTAF/..« <M, I) «*2* <Px*eZBOuN/ALAM) ♦TAUX)*ALAH 

1 aJ^tl-AMK A« A U / R.< N , I \*-' t B£T AF-^XAU4 

SN(N,I)=S(N*I) 

IE<hC0P.EU.MGO TO UQ 

TAljNs-tAU^-IAU Z#DZ — : 

TAuUsTAU 

CN(>.A)=C(1A)4TAU«0Z 

— ca< i a j =ct iax — _ — — ____________ — _ — __ — _ 

CZ(IA)sTAUN 

tau=taun 

GO — T o -L7 0 — 

I6f. TAu , '>=.5*(TAU04TAU 1 »TA|iZ*0Z| 

CZ ( £ A ) = T AUN 

CNxlA li.*S* {G 0 .U.A 4 iCilA).*,S* (XAU0 ^a-UN4-*0Z4 

l7r taN(.N«I)=S*JPT (S0F/(1 «*TAUN. o< *2) ) 

UN(N,I)=TAlN*wN{n»I) 

GO—TC -200 ; 

180 IF (LCOP.EO.I)GO TO IqO 
PN <N, I > sP (K * I ) 4P7*DZ 

SN ( N « li = 5 l K M ) ♦ SZ*U Z : 

UN(N,I)=U(N»I»*UZ*OZ 

V»N(N»I>=Wt\»I)*WZ«CZ 

GO— T©— 200 — — 
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Su p 02410 
-Su p 02420 
SUP02430 
Su p 02440 
-Su p 02450 
Su p 02460 
Su p 02470 
-Su p 02480 
SU p 02490 
Su p 0 2500 
-5u p 025 1 0 
SUP02520 
SUP02530 
-SUP 0 2 54 0 : 
Su p 02550 
SUP02560 
-SuP 02570 

Su p 0 2580 

SUP 02590 
-Sl)P02600 
Su p 0261 0 
SUP02620 
SUP02630 
SUP0264C 
Su p 02650 
-Su p 0266O 
Su p 02670 
•S|j p 02680 
-Su p 0 2690 
Su p 02700 
Su p 027l 0 
-Su p 02720 
Su p 0 2730 
Su p 02740 
Su p 02750 



19c — PN (N. I >=.5« (PO(N,I ) ♦,;(n*Ix*PZ*DZ>- 

sn(n« d=. 5* <so«N.n ♦« <N»n ♦sz^ozi 

UN(N t I) a.b* «UO(M*I ) ) ♦UZ*OZ) 

WN(IV t [ ) = .5« <WO(N, I>*t.r<N»I-y-*WZ*OZ>- 

200 CONTINUE 
21 0 CONTINUE 

8N f I )»HI| > 

BN (2 ) =CN ( 1 ) 

Z = ZN 

IF < 1 A .EQ.2T-CZ-M >=CZm 

8Z <2 ) =CZn 
00 250 1 = 1. IA 

NC=hA H4 +-1 

IF(*.EQ.1)nC=nC* 1 
8 ( I > s 8 N ( I ) 

C < I * sCN-( I > 

IF I^COP.E&.OGC TO 2->0 
• 00 220 Nsl.NC 

P-Olf* »-I T-sB 1NJ-I4 

UO (N, I ) =U (N » I ) 

W0(N,I)sw(N»I) 

22 O-SCUN * 11 s S (A *I-> 

23n PMl,1 )=pN(J,i ) 

W.N (1.1) sWN (3.1) 

SAU-L, U? SN (3.-U- 

UN (1.1 ) =-U\ (3,1) 

UN(2, 1 ) =0. 

— -00-24 0— N= 1 j NC - 

P (N . I ) =PN (N * I ) 

U (N * I > =UN (.\ * I ) 

W..( NJ I J SLW N. LlviXl-— 

S (N * T ) =SN ('• * I ) 

T<MI)=EXP(P{N,I)/GA*S(N.I)/GAMMA) 

2A-"-_R.lN » I ) = « ( I ) ♦X tN»4)_«4A( I T-.Rtm 

25 "■ GON'iiNUE 

IF (Lcop.eq.ijr'etuRn 

LOOP* T 

GO TC 7.' 

end 


SUP02760 

SUP02770 

SUP02780 

SUP02790 

SljP02800 

SUP02810 

SUP02820 

SuP02fi30 

SUP02840 

St>P02B50 

SUP02B60 

SUP02870 

SUP02880 

SUP02890 

SUP02900 

SUP02910 

SUP02920 

SUP02930 

StjP0?940 

SUP02950 
SUP02960 
— — SUP02970 
SUP02980 
SUP02990 
— — SHP03000 
SUP03010 
SuPO 3 620 

-SuPO 3 030 

SUP03OAO 
SuPO 3050 

SUP03060 

SUP01070 

SUP03080 

SUP03090 

SUP03100 
SUP031 10 

SUP03120 

SUP03130 

SUP03140 


-SU8H0UTlNE-OuTPUN(IBc<S,N7IRES.IT-yPc)- 


-O U n ^ 0 u I i HUN \ 1 7-1 * * * ' r * 

C0Pr i CN/BLM / NC(8) .PC (8 ) ,NC1 »NC2 • NCl .NC4 »NC5 .NC6 «NC7 .NC8 *PC1 .MC2 
l»MC.3.^C4,MCS,M*;6fMC7«MC8»NPEG(8) • NvC (2) »M«C ( 80.2) »NMAX »MpAX 
— 2 »6 AMm A » G A * GB.GC-, GO , Gr . GF , X (40.8) «-Y r l9.«) , XXP ( 130 ) 1 

3.YYP(13''.19).HM.XE.Yr.YA,xC,Ho.RD,nHFLO.TT.CC.EP,PTl 

4.SINTHE (20) *C0 SThE ( 2^ ) .9(20*19) M.SvMfLA.Ox (fl) ,DY (8) 

— COPP‘CN/HLK2 / P (l50»19 r *U( 150. 19) .V 050*19) ,S(150*1 9 ) *P I (15009) 

1 ,ui < i5() * 19) * VI ( 150 * lo) . SI ( 150. 19) »»S (8,2) .NF (8.2 > 

2 *MS( 8 , 2 ) ,MF (8.2) »TlMr«DT,K»J«XOTT(?)*OINF.OINFN*KOlVS,OTS 
-C0 M P‘CN/BL K 3/ YET ( 1 3c . i 9 ) . XCS ( 1 30) . Xp.X 1 f X2 , X3, YO » Y2 , ALP (4 ) »0D ( fc ) 

_ _ • , . . - - . 1 .... .nr ii ( 1 a 1 . U/*OQ ( 1 


0 
0 

— 0~ 

0 
0 
0~ 
0 
0 

-G- 

■ 1U U^v * » ' f A JVf | . fr'-Y nvf »w* , ' - -- - - - 

) .BE 1 ( 4 ) .LSLP.HU d 3o) „HL (130) .HC (1 35 ) »H(jPF ( 1 30 > »Hl-PB(130> »HCPR (130)0 
C0PMCN/HL^4/EP1 (20) .cP2(20> *EP3(105) .EP4I100) »EP5 (38) .EPt (38) 0 

4,NM ( 2(|) .NN2(20) ,Ml (=0) ,M2(20> .Ml C»8) »M4 (30^L'l (IOO) »L3(100) 0 

2.11 (lOo) .12(100) .JP(a0,19> *B(20) »8pR(2o) ° 

FPPES (PP) = E*P (PP) ,3, 0 


-V-^—UETU * V 

00020 
00030 
00040 
0oo50 
00060 
00.070 
00080 
00090 
00100 
00 110 
00120 
-0O130 
00140 
00150 




-FOE^S (PP.SS1 =ExP < (PP.SS>/GAMMA) 

FAS <PP«SS) =UF«SORT (EvP <PP/GA*SS/GAuMA> ) 
FACH (UU* W*AA) sSQRT <nll**2*VV*<*2) /Aa 
-FCP IPRE ) a2 . * ( PRE- 1 . ) ,<GAMMA*EM«o?>- 


FPTOT (P«E« iM ) =P«F* ( 1 ,*Gn*AM«*2 ) *«Ga 
1001 FORMAT <//4x «4hsTeP* Ic« 1 OX .SHTIMEsEi 2.4, 10x,3HOT«El2.4) 

1 002- FORMA T </4X. lOHREGlCN NO, I?) 

1 003 FORMAT </5X*-iHM=I3//9y*lHX»9XtlHY,q».lHP»9x,lHA*qX*lHUi<)X*lHVt9X 
1 *lHM,7Xf5rtSLOPE.GX,3uRH0.7X,2hCP,8»,lHS) 

1004-F0BMAT (I3« UF10.5» 

1005 F0rMaT(//4X*27hB0ATTaIL SURFACE FLnW FIELD /) 

1006 FORMAT (//4Xt2flHP0ATT A lL COMPLETE F ( OW FIELD /) 


_LOOI_F.ORMaT (/„/4X*27hNACEL;-E_COMPj.ETE-F-LAW-XIELO-/-) 

1Q08 FOrMaT(/4X,26hnac£LLc SURFACE FLOW FIELD //4X.5HSTEP 
l « I5»5X«5 hTIMEs,F6.3) 

_J.Q 09— FORMA T- (/4X • 1 CHOUTEfi — ftOWL— } — : 

1010 FORMAT(/4X»l0HINNEfi rOwL ) 

1011 F0RMaT(/4X.14HpLUME qOt/NOARV I 
_LQ L2— E-OfiM a T..I 1 3 X-t-lUX^UXt-luT-t-l-ixaHtUXl-X^-lMM-aOx-^H T A U^ lOX^ 2HCP*4-0X»2 hPT 0 

1/) 

1 0 1 3 FORMAT (4X»7F 12.4) 

_L0 14_:f. C K M a T-C / / 4 X >.13h G E 0 M E r RY- Tf S T A) — 

1 015 FORMAT (//4X»8HC0wL LtP / 1 3X . 1HX 1 1 1 x » lH Y* 1 nX , 3HTHE ♦ 1 OX * 1H0 ♦ 1 OX 
1 • 3H6PR ) 


— 0- 

-00160 

0 

00170 

0 

00180 

— 0- 

“OO 190 

0 

00200 

0 

00210 

— 0- 

00220 

0 

00230 

0 

00240 

— 0- 

01250 

0 

00260 

0 

00270 

— 0- 

-00280 

0 

00290 

0 

00300 

— 0- 

-003 10 

0 

0 0320 

0 

00330 

’-TO - 

00340 - 

0 

00350 


_U1 1 6_ F 0 R M A-T 1 1 3 X . 1 HX ♦ 1 1 X . LuY.lOX.3 HYPR+- 
1017 FORMAT (/4X,MHB0ATTAIi J 

IF (XTYPE.EC.3i GO TO 210 

IF-4LSYM. EQ..i4_G0--lC— ; *— : — 

IF(X.TYPE.EC.l) Go TO loo 
GO lO 4 

: .COMPLETE _F-LOW-.EJtELD-rtU.Tmi 

2 IF<lTYPE.EC.l)WRlTE( (t> »1005) 

IF ( 1 TYPE. EC. 2) WRITE (*.» 1006) 

.GO-t-O - 6 - 


0 

0 — 

0 

0 

0- 


00360 

00370 

00380 

00390 

00400 


WR.lJp(6» 100f) 

WRI-Tf (6. 1 0 0 A IK, TlME.iST 

-LRS = I.PEG __ 

IRF=IREG 


0 

00410 

0 

00420 

-0— 

00430- 

0 

00440 

0 

00450 

0— 

00460 - 

0 

00470 

0 

00480 

0 - 

00490 - 

0 

00500 

0 

0O510 

0— 

-00520 - 

0 

00530 


IF ( iREG.NE.9) 
IEs'"l 


GO TO iO 


0 

- 0 - 


O054O 

-0n550 


lO 


nO 


-iRF„=e ___ 

DO UQ LPEGsIRS.IRF 
IF(LrEG.EU.P) GO tc 

MCC=MC CLREC- 

NCC=NC (LREC) 

IF(LrEG.EO.P) GO tc iO 

IF <LREG.GT-.UAND.LREft.L-T.6l— OO—TP-PO- 

MSTAsl 

MFINsNTImES 

I F. < N T I ME S . E U. - M-MF I m=MC ( L R£G ) 

GO TC 4 (• 

20 MSTAsMCC^I-NtihES 

IF 1(NTImES.EU.-M-MSTaJs1— 

MFINaMCC 

IF(LsYM.EO.l) GO TC oO too 


-V -v V yvv 

0 

0 0570 

0 

00580 

D — 

00590 

0 

00600 

0 

00610 

-0— 

0«620 

0 

00630 

0 

0 0 64 0 

0- 

0O650 - 

0 

00660 

0 

00670 


0 

n 


00680 

00690 



GO— TC- AC— 

3n MSTial 

HFlNaNTlMES 

IF (.MTlMES.cU..i > -MFIm=M0-(L-HEG) 

IF {Lsya.nE . 1 ) GO TC iO 
IF (LsYM.EQ. 1 ) GO TC aO 

_.4C-.WRI If (6* 10C*> L«ES — — 1 

00 7 o m = mSTA ,mF IN 
*R I Tp ( 6 * 10: J)M 

JYY»7-(M.LREGl 

00 7(j N=1*NCC 
XX=*(N.LREO> 

nn=Nreg ilheG) *n — : : 

PRES=FPPES(P(NN,M) ) 

DEnSsFOENS (P(NN.M) «S(NN»M ) ) 

ASsiFaS (p (NN »M>.,S(NN,.*> 1 - : 

ACh=FACH(U(NN,N) ,V(Nm»M) .AS1 
ENT=S (NN«M) *1000 

CPsFcP (PRES). : 

IF (LcEG.NE.e) GO TC cO 

slcPe=o. 

UCaH.sU (NN»K).*C0STHELUi)--V(iJhLfAU«54.WiHE.tW4 

VCA h =U(NN»n)*SINTHE(v)*V(NN«M)*COStH£(N) 
IF(A0S(UCAR) .GT.l.E-oj SLOFE = VCAr/i .CAR 

XP s R ( N t N ) * C 0 S T H E ( N )_♦ 6 0 - RO 

YP=R (N,M) *SlNTHE (N) *uH 
GO To 60 

-50— CONTINUE ; 

XPsXxP (NN) -HD 
YP=YyP(NN»V) 

—SLOPE =0-. — - 

IF (Ars (U (NN *M) ) .gT. 1 E-R) $LOP£=V (N*'.M) /UInN.H) 

60 WRlTp (6,1 00*1 N'XP.YP,P«ES.AS.U(NN»w) *V(NN.M> ,ACH» 

_7q_cokTinue 

80 CONTINUE 

DUM = 1 ./ (GF-*tM*HH) 

ILthA-.EO. l )-PUN=2»*0uM/HH 

DO 9o L=1 i2 
NNM=NREG (3) 4 6 

— — ~I F- ( ); , E 0 • 2 ) - NN M s NR E G ( * T ♦ 10- 


SLOPEtDENS.CP.ENTO 


xem = x xp ( nnm ) : 

EMOOT = 0 . ’ 

EH I =o . 

. 00-85- M=P.hC3 ■- — ; 

EH I TsEhI 

OEL=VYP (nN^ .M) -YYP (NvM,m-1 ) 

DENS = F0ENS<P{NNH,M) »p<NNH,H>-1 

AS = F/is (P (NNM,m') »s (NN...M) ) 

ACH = FACH (U (NNM.M) ,v ImNm.M) .AS) 

EHl=5DENS*AS*AC h - — 

IF(La.EU.I) EMI=EMI«CYP (NNH.M) 

85 ENDOT=ENQOT* (EH I .Eh I ? ) «0EL/2 . 

EH0Ot=£M0OT*0UH ; 

WRltE (6.1010) xEH.EMnOT 

1018 FORMAT (//AX.2HXa,F7, 6 .5X.5HM00T=,Fa,A//) 133 



-9 ft— CO kT-I NUE 

J e 0 

RETURN 

-C SURFACE- FLCH-FIELD---NACEbL£— 

100 WRITE (61 10c«> K.TlME 
00 2 (jo L = 1»S 

GO- ic- < lio. 120.130 ♦ lio rlSOUtr 

1 1 0 WRI»E ( 6 » 1 OC 9 ) 

WRiJe (6 ,1012) 

LREGsS 

NSTA*1 

NF INsNCH/2* 1 

ms 1 : — — 

60 TO 160 
I2r LREG=7 

NST*=NNC4l-) 

NFlNsNC7 
Ms 1 

60 — IC— 160 

13(5 WPI]f (6,1010) 

WRITE (6, I0l2) 

LREGeft : 

NST£=NC8/2*1 

NFIN3NC8 

P=1 ; 

60 TO 160 
14o LREG=3 

NSTiJsNNC ( 2 ) 

NFINcNC3 

MsMC3 

GO-TC 160— - — 

150 WR I I e (6 .101 l) 

WRlTp (b.1012) 

LREGs6 

NST?=1 

NFIN=NC6 

16'; DO J. <3 6 NL = N5TA,MFlN 

n=nl 

1 F_(.L . E Q . 1 ) _NsNF I N-NL * 1— — 


0 01060 
0 01070 

-- 0 — OlOflO 
0 01090 

0 01100 
-0 — Oil 10 
0 0L120 

0 01130 

-O — 01140 
0 01150 

0 01160 
—0 — 01170 
0 01180 
0 01190 

-0 — 01200 
0 01210 
0 01220 
-O — 01230 
0 01240 

0 01250 

-0 01260 

0 01270 

0 01280 
— 0 — -0 1290- 
0 01300 

0 01310 

-0 — 01320 
0 01330 

0 01340 

-O' — 0 1350- 
0 01360 

0 01370 

-0 — 01380 - 
0 01390 

0 01400 

-0 01410- 

0 01420 

0 01430 

-Q — 01440 


.NNsNfiEG (LREG.L+N— 


U 


PRFSsFPRES (P(NN. m) 1 
ASsFaS (P (NN »M) ,5 (NN ,«) ) 

ACH = FACH (U (NN.M.) , V(Nm«M:)-,aS) 

CP = FCP <PRES> 

PTcTrFPToT (PPES.ACH) 

1 F-( L R E 6 . E U . 8 1- 60 -T O 0 

XPsXxP (NN) -RQ 
YPsYyP (NN* i '> 

SLOPE*:'.. 

IF(A0S(U(NN»M) ) ,gT . I ^ E- ft ) SLOPE* V (N n ,M ) /U ( mN » M ) 

GO Tc 180 

Uf._ XPs6 (N»M) «CUSTME (N) *n0-RD- 

YP=H (N.M) *SlNTHE (N) ♦qH 

slope=o. 

UCARsU (NN*PJj*C0STHE-L4X-V (^^HM»SlNiHE4M) 134 


0 01460 

0 01470 

-0—01480 
0 01490 

0 01500 

-0—01510 
0 01520 

0 01530 

-0 — 01540 
0 01550 

0 01560 

-O — 0157ft 
0 01580 

0 01590 

-0 01600 - 



VCaR*U(NN«v >*SINTH£ <m> *V(NN*M) «CriSTHE <N> 
IF<2eS(UCAc) .RT.l.E-oJ ?LOPE=VCAP/uCAq 

18 H - WR I f £ < 6 ♦ 1 0 1 J t~- *0, YP to«6SrAGrtrSL0P£vC^rP 

I9n CONTINUE 
200 continue 

.RETURN 

c geometry Test 
21 o WR I |e (6 * 1 0 1. 4 ) 

If..(Lc YM .£0 J— GO— 76— »25- 

WRITE (6*1015) 

00 F?0 Ns I * NCR 

I M£ ®PlI.°X.(N*8l *P.ll/2- ! 

» xpsB(N)«costhe ( N) ♦roIro 

Y P-e <N)»SJHTHE (N)*Hu 

220--WR I Te <6* 1 o t 3 >-*P^YP *tHE-*6 (TO-*^PR4N-'h 

225 00 *90 L*1 .4 

GO tc <230*240*250*2«,0) ,L 

234-1 F-( L S V-M. £ Q . i. )_G 0— T Q—i 33 

WRlTElejflO;:*) 

WR I Ie (6*1016) 

NST*sNNC ( U 

GO To 234 
232 WRlTp (6» 101 T) 

WR I Tf (6*100 9-)-— ^ : 

NST As 1 

234 NFlNsNCT 

! LR£^s7 

M 3 1 

GO To 27o 

2-4 0-1 E.(LS YM.EQ- 1 )-G0-T-C— >9o 

WRliE (6* .101.0) 
wpiTe (6* 101 P) 

flST5aNNC.t2.1. 

NFlN = NC3 

LREG=3 

M=M<:3 

. GO To 270 

25o IF (LSYM.EQ. 1) GO TC }9o 

-NS-T A s l— : 


NEIW=NC4 

L R E G a 4 
M*V^4 

- 60 - TC - 27 o * 

260 WRITE (6*101.1) 

WRiTe (6* 1 016) 

NST A» 1 ' 

NF I NaNC6 

LREif = 6 

270 DO 280 NsNSTA.NFiN 
NNsNRFG ( l R E u > *N 

-I F -( U t E Q • 1 •GR.L-. F0»4> -Mp O eHUP R< NN ) 

IF(L,EQ«2«CH,u.EQ*3) HpOsHLPR (NN) 
280 WPlTc (6* 1013) XXp (NN', *YYP (NN,M) »HPo 
-290- COnT INUE 

return 

End 135 


0 01610 
0 01620 

0- — 01630 

0 01640 

0 01650 

— 0 — 0)660 - 
0 01670 

0 01680 
— - 0—01690 
0 01700 

0 01710 

— 0—01720 
0 01730 

0 01740 

— 0 — 01750 
0 01760 

0 01770 

— 0 — 01780 - 
0 01790 

o oiaoo 
— 0 — 01810 - 
0 01820 
0 01830 

— 0 — 01840 - 
0 01850 

0 01860 
— 0 — 01070 - 
0 01880 
0 01890 

— 0—01900 
0 01910 

0 01920 

— 0- -01930 
0 01940 

0 01950 

— 0 — 01960 - 
0 01970 

0 01980 

— 0 — 01990 “ 
— 4 — 02000 
0 02010 
0 02020 
— 0 — 02030 
0 02040 

0 02050 

— 0—02060 
0 02070 

0 02080 

— 0 02090 

0 02100 
0 02110 
— 0 — 02120 
0 02130 

0 02140 

— 0 — 02150 - 
0 02160 
0 02170 



SUBROUTINE. SETNOS IJRt-TE)- 


COVf'CN/BLKl/NC (S) *HC <8) . NCI • NC2 » NC3 *NC4 « NC5 *NC6 t NC? »NC8 tMC \ *MC2 
1 .MC3fMC<nMCS,MUfe,MC7,HC«.NREG(8 ) ,NmC(2> «HmC(80»2) *NMAX»M*AX 
-2 » GAMMA »GA »G 8 »r;c »GO *Gp »GP « X <40 1 8)-. Y-f 1 9 « 8 ) » x X < l 3 0 > 


3»YV(l3!j»19) «HH.XE»YE'.TA»XC*RO.RD.R«PLO,T7 .CC»EM.PIi 
4»Si:MHE f20> *CO^THE (2* )iR(?0»iq)»LSvM«LA«DX(8)»DY(8) 
-COmM C n/UlK 4/EP1 (20) *pP2 (2n ) »EP3«1 0“ ) »EP4 n 00 > »EP5 08) ,EP6 (38) 


ltNA^ (20) *NN2(?0) »M1 (^0) »M2(20) *M3(i8) *M4(38) »L1 (100» »L3(100> 
2* II <100) »I2 < 1 "»o ) , JR Co 0.19) ♦B<20) .BdR(2o) 

— 1001 -FORMAT (4X.25HN0SE- RErt ION -INTERPOL A t 1 0N/4X, 3HB-C/8X f lHN , 2X 

1 *3HNM «2X*3HNN^ «3X»2uMi ,3X*2HM2t 1 X . 4h I 5 1 T ,ftX « 3HEP1 « 9X « 3HEP2 ) 

1 002 FORMAT (40flHlNjTERlCR/8x, 1 hN*4Xi1hM, 3X.2HI) t 3X * 2R 1 2 , 3X » 2HL 1 
— K3X*2HL3,6X ,3 heP3»9XOhFP4 ) 

1003 FORMAT (4X*4MX = R 0/8X*,HM* 3X *2HmR, 2X.3HNNC *4 X*lHNt3X,2HM3»3X 
1 *2HM4«6X*3htP5»qX«3HcP6) 

14)04 — F.ORMA T-(..4.X« 6 L&. pE 12^4', 

c set up for interpolation at r=c 
CC sRo 

00 — I(j— Ns 1 -*NC 8 


THE»PII*X <M8) *PlI/2 t 
CALI; WALL(8*THE«BB.8nPR»3) 
-fllN La 89 ' 


BPR ( N ) sB0PB 
COSTHE(N)=COS(THE) 

-SlNlHE CN1_=5 INiTHEJ - , 

DO io M=l • C 8 
R(N«R)=Y (M,tjj* (cc-BBi *BB 

-ICl CONT inue _____ 

NA8=NC8-1 
NMs;'- C 8/2t1 . 

IELX>)RI T£*NE • Q j W R I_I E. * 6 .. I o Ql L 

00 *00 N= 2 iNar 
IF (N.EQ.NM) GO TO 13? 
-IHE=PII*X(N*8)*PH/2- — 


20 


XA=CC*COS(TME)*Ro 
YAsCc*SIN(THE) tw H . 
-IF-(Th E«-G£.£ l£l_GO TO. 2q. 
LRE<J = 7 

MRE 6 s 1 .... 

—GO — T.C — 3_ ) 


LREf =3 
mreGs 2 

lO.NC.C =NC(LREC) 

00 5? L=l*NCc 
LLsNREG <LPEB) * l 
LF-l * x (LU_.e r.XA L 


So con' INUE 
6 C L11=L-1 

XBaAx (LL1) 


_G£ t Q . 60 - 


RB1= (XB 
_A 1 . 1 _t.XB_ 


■RO)/Cns (THE) 


IF (MREG.EOiU GO TO pO 

00 JO .Ms I i”C7 

_RRL= — SQRTX ( All - R 0) •*>*_ ( YY (LL1 -* H ) -H^) ««2) 


-SElOOOlO 
SETO 0020 
SET00030 
-SEI00040 
SET00050 

set oo o 6 o 
SET0O070 
SE 100080 
SeTO.ooRO 
-SeTooIOO 
SeTOoUO 
SET00120 
-SET0O130 
SET00140 
SET0O150 
SE TO H60 
SET00170 
SET0H80 
-5ET0H90 
SET00200 
SET00210 
— SE^0o220 
SET01230 
SEl 0 1240 
— SET 01250 
SET 00260 
SET00270 
-SE701280 

set oo2*?o 
se to o 300 

-SET01310 
SET 00320 
SeT 00330 
-SET01340 
SET01350 
SET01360 
-SE-T 01370 
SET00380 
SeT 0 1 390 


-SET01400 
SET01410 
SE TO ^ 420 
-SE TO 0430 
SET01440 
SeT 0 1450 

-SE To 1460 
SE TO 14 70 
SET01480 
-Set 0 o 490 
SfToisoo 
SET00510 

-Se to i52o 
Set 01530 
Set 01540 

-SET01550 



I F. I P RJL.G T. » C „? L jGQ __T£_i5. 

70 CONTINUE 
75 mm=M 


GO Tc - 90 : 

00 00 *4 LMai* M C3 
M=vC3-LM*l 

RR ) .SORT U * 1 1 -RO ) * *? ♦ < Y Y ( LL 1 «*2 > 

IF (R<?1 ,GT»CC) GO TO «5 
84 CONtlNUE 

fiS_MMsM : 

9n THE A laATAN2 ( YY (LLl «'Mm) -HH.X1 1-RO) 

I F < T H E 1 l.LT.C.) THE17 a2.*PII*THEli 

— LL2=LLUl — 

IF (THE 1 1 .LT.ThE.aNC.mREG.EO.1) LL2 = LL1-1 
IF ( The 1 I .GT .ThE.aNC.uREG.EG.2) LL2-LL1-1 

IF_C L l2 . LT . N-HE G (LREG ) *T4_-LL2=LUU4 

X22 = X X <LL2) 

Y22=YY (LL2* m M) 

Y 4 1 - Y Y-< LL 1 • Mm ) 

0EL=(Y22-YU)/(X22-X?I) 

RB 2 S <Y 1 l-Hh*UEL*(Rn- Y ll) ) / (SIN (T hE ) -0EL*COS (THE) ) 

1 E-( L 1-1 .£0. l)~RR2=lC-» 

IF (R82.GE.CC) G o TC nl 
MMsfM- (2*MRtG-3) 

G0-Tc_9. 

9 l CONTINUE 
I S I T= 1 

IE_C K 82 »LT • F.W-14 — I-S-Ila? 

IF(islT.E'J.2) GO TC 720 
IF (MCEG.EQ.2) GO TC T 05 

00 -1 6 0 -Ms i . «C 1 — 

IF (Yy (LL l *Y> .GT.YA) „0 TO 1C4 
100 CON • INUE 

1.04. MMl=M . ^ 

MM2=M-1 


GO To 118 

10S-00-1) n-LMSU..»MC:7— : 

M=mC7-l^*1 

IF(Yy(LL1* v ).LT.YA) rsO TO 115 

ULO-COn i INUE 

115 MM 1 = M 

MM2 = M* i 

lXa_Y.Bsf1H.RBlF.Sl NIT HE-)— 1 


EP1 <N ) = ( yB-TY (LL l iMMT ) ) / ( YY (LL l »MM?) -YY (LLl »MMl ) ) 
NN1 (N) =LLl 

NN2(A) =LLl : 

Ml (N) a MM 1 
M2 (N)=Mm2 

R B=RB1 

GO To 150 

120 XB=R82°C0S (THE) .RO 

EPl-U) stxri-'Xll )/ <X22.X114 : 

NN\<M=LLl 

NN2(n1=LL2 

Ml-C)sMM ^ ™ 


__.SeT0*»560 

SETOT570 
SeT 0">580 

— set 01590 
SFT01600 
SFT01610 

SFT01620 

SFT01630 
SET 01640 

SET01650 

SET01660 
SeT 01670 
— SET01680 
SeT 01690 
SET00700 

SET0O710- 

SfT 01720 
SET01730 

SET01740 

SET00750 
SeT 01760 

SET00770 

SET 01780 
SET01790 
— SET01800 
SET01810 
SET01820 

SE-T01830 

SET 01040 
SET00850 
— set oi86o 
SET 0 1870 
SET 01880 

S E-T 01890 

SeT 01900 
SET0O910 

SET 01920 

SeT 01930 
SET 01940 

-SET 00950 

S E T 0 1 9 6 0 
SETO 0970 

SET0O980 

SET 0 0990 
SET01000 

— set o i o to 

SeTOI 020 
SeTO 1030 

SET 01 04 0 

SET01050 

SET01060 

SFTO l 0 70 

SeTOIOBO 
SF T 0 1 090 
— SeTOI 100 



J424NV3MM 


R8*He2 
60 Tc 150 

+30-- N N 1< M =NREG <7 »*} 
NN2 (M =NREG ( 7 » ♦! 


Ml (N) =1 

M2 (Plj.32 

EP1IM»C. 

RB=Cc * 

l-S/v_RC»~tcC-B (N) ) *Y(Mc8~}v8^*&^) 

EP2<M = <CC-«C>/<RB-Rr> 

IF(JfilTE.NE.O) 

fWR l lE <6 HCOM— H^NN-1 (MKNN^?+N4^M^HNi-rM2-fWl-»-!-SIT-,ERHN->-»EP2 (W+ 

200 CONTINUE 

C SET UP EOR INTERIOR oOlNT INTERPOLATION 

jj=o : — : 

IF (JRITE.NE*0) WRITE 1(6*1002) 

00 *00 L I = 1 * 2 

IF-lL I .EO.^U-GO-TO-2 - ^ 


LRe|J = 7 
GO 10 220 

2Lc_LREG=3 : 

220 NCC=NC(LREG) 

NCC=MC(LREC> 

0O-3S0-N=2 »NCC— ^ 

NR = < L I - 1 > *N*AX*.N 
NN=NrEG (L«FG) *N 

IE< * X<NN>.G I-.R o l—GC-.; 0-360 

00 ^40 LM = 1 *MCC 
MaLM 

IF-tLUEQ.2 1 — M=MCC-®'W^-*-lr— — 

XAs a X(NN) 

YA=Ty (NN,M) 

R AsScRT l (AA.-JR04 iA-HH4-?-*2J 

IF(RA.GT.CC> GO TO 3a5 
THE«=ATAN2 ( YA-HH,XA-oO) 

IF-<THEA-.LT-«o»4-TWEA—a?«*PII-*^HEA 

JJaJj+l 

JR(Nr,M) s JJ 

00— f 5 0— L= I-t.NCR 

THE S P 1 1 * X <L»8) ♦PlI/2. 

IF (ThE.GT.TMea', GO Tn 260 

25_0_jCQN-I-I NUE- „ 

26o L3<Jj)=L*NcEG(0) 

LI < jj> =L*NRtG(a)-l 

IHE-3=THt 7 

THE a = PI I*X < I — 1 ♦B>+PIt/2, 

DO 270 I=1»MC0 


RRa lCC-BlL-D ) * Y < J-* fl * ♦ B-( L- 1 > — 

IF (RR.GT.Rfl ) GO TO 2a 0 
27 n CONTINUE 

28q I F_l J . E Q . 1 J._ 1=2 

12 ( Jj)=I 
II (Jj 1=1-1 

RR2--tCC«B-l.L>l-M LI * o L+83L---U- 
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- - setouio 

SeTOI 120 
SETOl 130 

SET0M40 

SeTOI 150 

SeTOI 160 

SETOI 170 

5ET01180 
SeTOI 190 
- — SeTOI 200 
5ET01210 
SET01220 
— SE TO 1230 
SeTO 1240 
SeTOI 250 

SETOI 260 

SeTOI 270 
SETOI 2fl0 
—SETOI 290 
SeTOI 300 
SeTO 1 3 1 0 

SE-T-0 1 320 

SETOI 330 
SET 01340 
SETOI 350 
SeTOI 360 
SET01370 

SETOI 380 

SETOI 390 
SET01400 
— SET01410 
SETO 1420 
SETO 1430 

SET 01440 

SETO 1450 
SETO 1460 
—SeTOI 4 70 
SE.T 01480 
SET01490 
— SE-T 01500 
SET01510 
SET 01520 
—SeTO (530 
SET 01540 
SET01550 
— SETO 1560 
SET01570 
SeTO 1580 
— SeTOI 590 

SeTO 1600 
SETO l 6 1 0 
— SE-T 01620 
SETOI 630 
SET01640 
— SET01650 


RR1 = ( CC- 0 1.L'ln * Y (I -T»_8 L* 8 (L-J.) 

RR3 S (CC-B<L> )*y(T-1«o)*B(L) 

RB4* (CC-fl (L > > *Y < I * a ) *0 (L ) 

EP -3 tjj) = (ThLA-THEl ) /-',.ThF3— THE-14 — 

EP4 ( JJ> -.5* ( (WA-RR2 ) / (RP1-KP2) ♦ (PA.RR4 ) / (PR 3 .RR 4 ) ) 

IF< JRITE.NF.O) WRITE >6,100«)N,M,I1 f JJ) ,I2(JJ) #L1 (Jj) »L3 <JJ> 

1 * E D 3 ( J J ) .EP4_(jj) — i 

34o CONTINUE 
345 MMC <NP«LI>=* 

35i COnT INUE. , 

360 NNC<lI)=N 
400 CONTINUE 

jC SEtUP^EOR_.l-NTERPOLAiIONLALQNG_X=R* 

IF ( JRITE.NE *0) WRITE ,6,1003) 

DO .666 LI»1 *2 

NRC=NNClLl-l 

IF (tI.EQ.2) GO TO 51* 

N=1 

N-l =2 

LREG=.7 
GO To 520 

— 5Lo— n=n ^6 : 

MsNc«-I 

LREG=3 

52-f- -MCC = MC (LREG-T - : 

NNRCsNREG (i*eg) ♦nrc 
XAA=RO 

00-59 0 MsLt.MCfl - - 

MR= ( l I - l I* vM Ak*« 

YAA = R (N*M) oSIniTHE (N ) *HH 

— — Y CC = R < N 1 , M ) * S I N T H E ( N r ) *MH 

XCC = R(Nl,M)*CO«;THE(NT>*RO 
X 8 R=XX(NnRC) 

— .- Y 8 B= UXUB-XCC) /-(XAA-vCCI 4 «UJYjAA«-Y-CCi-*TTCG- ■ — 

IF ( I * 1 • E 0 , 1 ) GO TO 550 
00 530 MM= 1 «MCC 

IF ( Y y ( NNRC«MM) t GT. YBo)_00— 10—540 

53o C0N]INUE 
54o MM 4 =rm 

MM 3 = K M- 1 

Y44=YY (NNWC*MM4) 

Y33=YY (NNWC»Mm3) 

£P.5(M.R) = ( YEP-Y33) /.(Ya.4-733-1 - 

EP 6 (MR) = (XAA-xcC) / (XqB-XCC) 

M3 ( Mr ) =Mm4 

M4.C.MR) =MM3 : 

IF (t'RITE.NE .0) 

IwRIJe (6 * 1 Ot, 4 ) MtMRfNoC, N,m3 (MR) ,m4 jMH) ,EP5 (MR) ,£P 6 (MR) 

GQj.T.C.-59.0 : 

55o IF ( L I , E<J ,2 ) GO TO 56i 
M3 (Mr ) = 1 

M4.(Mr ) =2; 

GO ]C 570 
56o M3 (Mr)=mc3 

M4( MR.) = M C 3 -J ; ; 
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-SETOI66O 

SET01670 
SET016BO 
— SE TO 1 690 
SeTO 1700 
SeTo 1 710 
-SET01720 
SeTo 1 73n 
SE^O 1740 
— SE To l 750 
SeTo 1760 
SeTo 1770 

-SEl0l7«0 

SeTo 1 790 
SET01800 

-SeTo 1010 
SeTo 1 b2o 
SeTo 1830 
SE-T01840 
SeTo 1850 
SET01860 
-Set oi87o 
SeTo 1880 
SeTo 1 890 
-SeTo i9oo 
SET01910 

set 01 920 
—SeTo 1-930 

SET 0 1 940 
SETC1950 
-SeTo i960 
SETO 1970 
SeTo 1980 
seto 1990 

SET 02000 
SET02010 
-SET02020 
SET02030 
SeT 02040 
-Set 02 n5o 
SET02060 
SET02070 

set 02 080 

SETO 2 090 
SET02100 
-SET021 10 
SET02120 
SET02130 
-SeT 02140 
SET 02150 
SET02160 
-Sp t 02 1 70 
SET02180 
SET02190 
-SET02200 



5 7r FP fttMRl = o • 

EP6 (*R) = (XflA-xCCJ/ (XqB-XCC) 
590 CONTINUE 

-&DC—CONI-I NUE — 


RETURN 

ENO 


_SU HR CU T I N E S.T RECHtM R e.G )_ 


COmMcN/BLK^NC (R) ♦ MC f8) *NCl •NC2«NC3»NC4*NC5*NC6.NC7*NCfi*MC1 *MC2 
1 »MC3 *MC4 ,MCP*MC6,MC7 ,Mcfl»Nfi£G (8) ,NmC (2)»MmC<80*2) *NMAX ,MM Ax 
-Z* G AMM a « G A « GO « GC ♦ GO . Gtr • GF ♦ X-( 4<) ♦ 8 ) « Y-f 1 9 ♦ ft > ♦ X XP fl 3 qT- 


3*YYP (13c»l9> fHHtXE* Yc*YA,XC*Ro»RO*cMFLO«TT*CC*EM*PlI 
4,SINTHE (20) *COSTHE (2*) ,R(2O*19 )*lSvM«LA*0X(8)*DY(8) 
_COmMXN/BlK2/-YETL< i3om*> *XCS41301.X4*X1 .X2,X3«Y0»72»ALP(4I *00(61- 

I r» ^ I i A i - I Cl U . . < « *3 .. I 1 Li. i « « I I ./» i » •>* I . . .. .nn . ■ . ^ . . . r»_ . < 1 . > > . . • . 


wv— ■ i -'V — T ''F T "V w»V r T r. ’TT "IF « *. » ft JF T «l. r 1^7 1 C 1 

1 »HE> (4) *lSLP«NU <13o) .Ml ( 130) *HC <1 3') *HUPR (130) *HLPR(130) »HCPR (130) 
OlPtNSION AL. ( 2 ) * ERR ( 9 ) 

0 AT* -LAMP / - /. 


1001 E0RMAT(/4X,2lHSTRECHfNG PARAMETERS /4x,4w X0=*El2.4*3X 

l * 4H Xl = »El2»4,3x,4ft v 2 = ,E 1 2 , 4 , 3 X . 4 U X3 = « E 1 2 . 4 . 3x * 4 m Y0=*E12.4 
_2-Ox»AH_Y2=itl2-.4/4x,ZHALl = »El2«4«3v-,4HAL2a,E12.4-»3xt4HAL3 a - 


3tEl2.4,3x»4NAL4=,El2.4/4x,4HBEl=,El2.4,3X,4HBE2=*El2.4,3X 

4»4HPE3=»E12.4,3X,4MBr4 =f El2.4) 

-E-l (F:*AL «a£-)-=.5« (4. * TftNH (AL*TF--.5)-VveE)- 


F2<G,AL«RE)=2.*8E/ (Ai *(V.-(2.*G-1 . )«*2*BE**2) ) 
G1 (F «A1 *A2) B (Ai *A2«Ai 0G(f > )»AL0G(F» 

-G2 < F-» A 1 . A 2) =4 A 1 2 . * A p « A L 0& < F-l 1 /F 

IF (Lamp. nE. oi go to T02 ‘ 
stpEtching Parameter*: 

-OO — 4- L=1 * 

IF (L.EQ.2) 

0EL=0X ( 1 ) 

-01=00(3)-— 


«0 TO 1 


02=00 (4) 

GO To 2 
O_0EL=0*(U— 
01=00(5) 
02 = 00 ( 6 ) 
-2-OElUI.^OEL 


DLOGsALOG (C tL) 
01lOG=ALOG (OEL1) 
OE-T ^CLOG-0 1L0G- 


, DUMls <02*01 LOG/OLOG-r. 1«OLOG/01 LOG) /OET 
DUm 2= (01/01 LOG-02/0LnG) /OET 

_J F_( L . E Q . 2 )_0 0 _ TO 3 

XOsOLMl 

X2=0lm2 

-GO TC 4__ : 

yo=pumi 
Y 2 = 6 U M 2 

-COnT-INUE : ! 1 

X1 = X0 
X3*X2 

.00 S-Ls 1*4 ■ ■ 

ALP ( L ) 3 C • 

8ET«L)=0. 

-00 iO_L=l-..4 
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-SE-T022 1 0 
SET02220 
SE T 02230 
-SET02240 
SET02250 
SET 02260 

-STROnolO 
StPOoo 20 
STR00030 
-ST«0n040 
STR00050 
STRO 0060 
STR04070 
STR00080 
STR00090 
— S TR 0 0 1 0 0 
STR00110 
STR00120 
- STR00130 
STR0O14O 
STR00150 
-STROM 60 
STR00170 
STR00180 
-STR0O190 
ST«00200 
STR00210 
-STR00220 
STR00230 
STR0O240 
STR00250 
STR00260 
STR00270 
STR00280 
STR00290 
STR00300 
-STR0O31O 
STR00320 
STR00330 
STR00340- 
STRO 0350 
STR00360 
-STR0O37O 
STR00380 
STR0^390 
-STR00400 
STR0O410 
STR00420 
-STR00430 
STR0O440 
STR00450 
-STR00460 
STR00470 
STR0O480 
STR00490 



GO To ( 11* 12* 13* H)',L STRO 0500 

11 IF(LSYM.EU.l) GO TC 20 STR0O510 

DEl = 0Y(2) ■ STR0O520 

OIm = 0EL*HH StRO 0530 

IF (00 (1 ) .6T . ,q«OlM) nO(l)3.9*DIM STR0O540 

0UM = 2.*D0 (-l-i /HH-l • STA00550 


60— T 0 — 1 5 : 

12 IF(LSYM.EQ.l) GO TC 20 
DEl=0X<3) 

01 M = 0 E L * X C— r- 

IF (Uo (2>-*GT «: .9*0 IM) pO (?) = .9*0 IM 
OUM=2.*00 <21 /xC-i , 

GO -TO -15 — 

13 IF (LSYM.EQ.l ) GO TC 2o 

DEl=OX (4) 

OIuscEL* <*E^XC) 

IF <Uc (2) *GT • ,9*0 IM) pO (?) a.9<*0IM 
0UM = 2.*D0 <2 ) / ( vE-XC ) -1 . 

GO -T O-— 1 5 

14 DEl 3 0X(7> 

OIm=oEL*X£ 

IF (Do < 2 ) .GT • .Q«r>TM>— nO-(3-)-»»-9*QIM 

DUM=p.*0D<2 > /XE-1 * 

15 A L ( 1 ) = • 5 

AL-<2>a.6 : — 

MEal 

KIP=1 

-IG-EFsTANH (AL(ME)fl <oEL«-5)4VTANH(-.5«Af-4>‘E-H 

ERR (PE)=EF-UUM 

IF(Ars<ERR(ME)), l T. 1.E-3) GO TO la 

1 F-< Mp . E u . 2 TO 0 -T 0—4 7 — 

ME*2 

GO To 16 

— 17 — ALR”r = A U ( .1.1 “E RR.(-l-) •« ( a L 4 2 T-AL4444V4c rr ^ ERR (-1-14 
ALRAr=ABS<ALbaR) 

AL(1)=Al<2) 

a L(2 >=alha p_ : 

ERR 1 1 ) =ER« (O 
KIP=KIP*1 

IF_(KI P .LE *2 0 1GC_T0_J?6 : 

CALL EXIT ‘ 

18 ALP <t ) =AL ( v t ) 

BET (L ) =TANK ( *S«AlP-<C+4- 

20 CONTINUE 
LAmPs 1 

— - — WRITE (6* lOl'l )X0 »Xl »Xa* X3 • YJ * Y2 t Ai.P-.BE T- 

STRETCHING DEPlVATlVrS AND BODY GErMETpY 
102 IF (MrEG.NE.«) GO TC 704 

MS T As 1 ‘ 

MF I Ns 7 
GO To 1 00 

L04_MST”sMPEG 

MFINaMHEG 

108 DO 2po LREC=m<;T A ,MF Im 

-NCC = NC ( LHEG ) : 141 ’ 


STR00560 

STR00570 
STRO 0580 

STRO 0590 

STRO 0600 
STR00610 

STRO 0620 

STR00630 

STR0O640 

STR0O650 

STRO 0660 
STR00670 

STR00680 

STR0O690 

STR0O700 

-STR00710 

STR00720 

STR00730 

S T R 0 0 7 4 0 

STR0O750 
STR0O760 
— STR0O770 
STR0O780 
STR00790 

STR0O800 

STRO 0810 
STR00820 

STR00830 

STRO 0840 
ST R 00850 
— - STR 0 O 860 - 
STRO 0870 
STRO 0880 
STR00890 

STRO ogoo 
St r 00910 
— STR0O920 
STR0O930 
ST R Oo940 

STRO )950 

STRO 0960 
STR0O970 

STR0O980 

STR00990 

STROlOOO 

STROlOlO 

STR01620 

STR01030 

stroioao 



MCC s PC(LRtG» 

00 1<J5 Nal.NCC 
NNaNfiEGtLR^) 

XXa* {N*LPE0> 

00 1<J5 Mal«MCC 
JY.Y.aXi_M ±L 8 E CJ 


STR01650 

STR01060 

STR01070 

STROI 08 O 

STR01090 

stroiioo 


M_X0._(110»i20 tl3Qj_U-0-«-LSiLti-&Q tl70) ,LRE& 

1 1 0 IF(N.EQ.UGO TO 111 
XXP(NN)=Gl <AX«xO.X2) 

XCS-< NN) = L «./<» 2-C XX *XQ »-»2l— 

GO Tc H2 

111 XXP 1nN)=-5c£>, 

XCS l\N) 

112 IF (T'.EO.MCClGO TO 111 
YYPUNfM) ayft-Gl (1 .-YCtY0*Y2) 

YE-I-i N N « M L3---U /G2 1 1 • -iY.T0.X2J 

GO Tc 1U 

113 YYP<NN*M)=500, 

YET< «M.= , * : 

114 C A(_h WALLU »XXP(NN) 

HU(NM=AH 

HUPfi < NN.).* AfeJiR 

GO TO 190 

I2p IF(LSYM.EO.l) GO TC ?90 

IF-(N .EG . 1 ) GO -Tq 121 

XXP (NN) =Gl (XX.X0.X2) 
XCS<NN)=l./G2(xX.X0fv2) 

_ _ GO Tc 1 22 - 

121 XXP 

XCS1nN)=0. ‘ 

122_E4A = F1 <YY*ALP<14*8ET4tM 

YYP (NN.M =YA«eta 
YET(KN*M)sF2(e;Ta*ALPH) »BET (1) ) 

____.CALl--ilALL(Z«JCXP(NNl f aH,AKPR»-U 

hl«nm=ah 

HL P« (NN) =A)-PR 

CALL -WALL < 2*X XPTNN )-^i H,^HPR#2-) 

HC(Nn)=AH 
HCPR (NN) sAhPR 

GO -T 0—19 0 

130 IF (LSYM.EG. 1) 60 TC ?90 
CSI*F-1 (XX»ALP(2) *8ET'<2) ) . 

XXP (fvN) sXC«CSI— — 

XCS InN)*F2(.CSI.AlP(2i *RET(2) ) 

CALh WALL<3 *XxP(NN) »aH,ahPR*1) 

HL (Nn ) aAH — : L 

HLPR (NN) =AuPR 

C-ALh WALL <3 *XXP (NN) »XH,AHPR*2) 

HC (NN ) =AH 

HCPR (NN) aALPR 

ETA=F1 (YYtALP(i) ,8ET' ( ln 

_YYP (NNtM) =ETA*HL (NN) 

YET <NN.ML=F2(£tA,ALP>1) .BETID ) 

GO To 190 

1A0_I E ( L STH . EQ . Ll _G0_XC_‘i-9 O 142 


— STROIUO- 
STR0U20 
STR01130 
— STRO 1 140 
STR01150 
STR0U60 
—STRO 1170 
STR01180 
STR01 190 
— STR01200- 
STR01210 
STR01220 
— - STR01230 
STRO 1240 
STRO 1250 
— S-TR.0'1260- 
STR01270 
STR01280 
—STR01290- 
STRO 1 300 
STR01310 
— STR01320 
STR01330 
STR01340 
—STRO 1350 
STR01360 
STRO 1370 
— STR01380 
STR01390 
STRO 1400 
—STRO 14 10 
STR01420 
- STR01430 
— STR01440 
STRO 1450 
STR01460 
— ST«01470 
STRO 1480 
STR01490 
— STR01500 
STRO 1510 
STRO 1520 
— STR01530 
STRO 1540 
STRO 1550 
—STRO 1560 
STR01570 
STR01580 
— STR01590 



CS I B F 1 ( X X ♦ A L p ( 3 ) ,B£T',3) ) 
XXP<NN)=XC* <XF*Xc>«CcI 
J(CS (NN) =F2 (CSI ,A(_P .t3V*BEX.l3).l_ 
CALh WALL H*XXP (NN) m H, AHPfi»l) 
HL (NN) =Ah 

.HLPRTNN )=AhPR 


ST«01600 
STR01610 
-STR01620 
S T R 0 1 6 3 0 
STR016A0 
-ST r 01650 


— CALi. WALL (AtXXP (NN) ,aH,AHP(U2) 

HC (NN) sAH 
HCPR (NN) sAhPR 

ETA=F1 (YY.ALPd ) ,BET,1>> 

YYP (NNfM) = HC (NN) ♦ (HL (NN) -HC (NN) ) •E'tA 
YET (NNtM) =F«!(EtA*ALP/ H»B£T <11 ) 

GO- To > 90 

15o IF (LSYM.EO.l ) GO TC ? R 0 
IF(N,EQ.NCC>GO to 1 5 T 

—XXP (NN) =x£-0 1 (1,-XX, yl t X3> 

XCS <NN) = -l ./G2 ( 1 .-XX', Xl ,X3) 

GO To 152 

1 5-i— x xp Inn > = 5 o c- ♦ - 

xcs Inn) s o • 

1S2 CALL WALL (5 »XXP (NN) »SH f AHPR«l) 

HL-(NN )-=AH 

HLPR ( NN ) = At- PR 

CALL WALL (5* XXP (NN) ,aH,AHPR,2) 

HC < NN )-= AH 

HCPR (NN) = AHPR 

ETA = Fl (YYtALP(l) ,BET',1) ) 

Y-YP I N N * M ) = h C .( Nn )- *-( -Y F «HC-( NN-14 -*E-T-A 

YET‘NN.N)=F^(ETA,ALP,1) ,HET (1) ) 

GO To 190 

16 C._I F ( N . E 0 . NC C ) G 0. -tO_l 67 

XXP(NN)*'XE-G1 (l.-XX*yl,X3) 

XCS<NN) =-l ./G2 ( 1 .-XX'.Xl ,X3) 

__GO-Tn.-l 62 L. : 

161 XXP TnN.) =500 • 

XCS<NN)=0. 

L62—CALL -PLUBO ( XXP_tNNU Au t AHP.BJ 

HU(NN)=AH 
HUPR (NN) =AwPR 

IF.( M, EQ . MCC ).GO — To_l6^ 

YYP(NN*M)=WU(NN)-G1 <7.-YY»Yo»y2) 

YET «NN*M) =-l ,/r.2 (1 .-vY,yo«Y2) 

GO -T 0-190 : 

163 YYP <NNiM) =500. 

YET (NN *M ) =0 • 

GO -To 190 

1 7 0 CS I *F T ( XX * fl LP ( A ) * BE T < 4 ) ) 

XXP (NN) =Xt«CSI 

XCS<NN)=F2(.<-SI,AlP-(A'> .BETtAJ-J 

CALL WALL (7 * XXP (nN) t aH, AHPRi l ) 

HU (NN ) =AH 

HUPR ( NN) =AhPR 

IF (N.EQ.MCC)GO TO 177 

YYP (NN»M)=HU(NN)-Gl ( 7 . -YY . Y 0 * Y2 ) 

Y E_T INN i M)=»l. /-62 (1 v Y» yO « t2 143 


~ -STR01660 
STRO 1670 
STRO 1600 
— STR01690 
STRO 1700 
STR01710 
— STR01720- 
StRO 1 730 
STR017A0 
— STR01750 
STRO 1760 
STRO 1770 
— STRO 1780 
STR01790 
SyROlBOO 
— STR01810 
STR01820 
STR01B30 
— STR01BA0 
STR01850 
STR01B60 
— STR01870- 
STR01880 
STR01890 
—STRO 1900 
STR01910 
STRO 1 920 
— STRO 1 930 
STR019A0 
STR01950 
— STR01960 
STR01970 
STR01980 
— STR01990 
STR02000 
STR02010 
— STR02020 
STR02030 
STR02040 
— STR02050 
STR02060 
STR02070 
— STR02080 
STR02090 
STR02100 

-STR021 10 
STR02120 
STR02130 
— STR02140 



GO To 190 
1 7 1 YYP <NN»h) =500 

YET UM«M)_=:i , 

19^ CONTINUE , 

195 CONTINUE 

— 2Ho cqn.t_i.nue ; 


STR02150 

ST«02160 

STR02170- 

STR02180 

STR02190 

STR02200 


HExUftAL_ : , STR02210- 

END STR02220 


SUBROUTINE WALL (lREG',X,Y,yP*LP» WAL0O010 

COM'’ 1 CN/bLKl / NC (8) »PC ftt) tNCl »NC2*NC1»NCA»NC5»NC6«NC7*NC8*PC1 »MC2 WAL07020 
l »MC3,MC4,MCb,MC6,KC7.McB*NREG(8) ,NvC (2) ,MmC (80 *2 > »NMAX,MpAx waL 0*030 

„2«r,AMpA.GA tGB«GC.GO,Ge- tGF-.XA(4o»8V ♦ v Y ( 1 9 *8 1 -. X XP < 130) — — WaL0OQ40 

3* YYP (13C * 17) »HHiXF« Yr» Y 4 , xC » «7 • PI7 , oMFLO « Tt , CC * E* * P 1 1 WAL07Q50 

4«SINTHE (20) »COSThE (2i J ,P(?0tl9) *1. SvM» LA,0» (8) »OY (8) WaL0*06O 

—DIMENSION 1«E (9i ,BN-0ct9 )-*bPN0SI9)-*AMl (9>-* AN2 ( 91-^X0l>T49MY0uT (91-* — WAL 07070 
1YP0UT (9) « LF UNO ( 9 ) tAO? (9) ,A,02<9) * X Iv ( 9 ) , Y In (9 ) , YP I N (9 ) »LFLNI (9) WALO 0080 
2«All (9) »AI2 <9) ,A(20) WAL0O090 


— ACUP4-XL-* XR » Y L ♦ YR *-Y-L P4TRP-) =-I4-YtP-*-Y-Hpl x R )»2 c* YR ) > / (XL-XR) WaLOOIOO 

1 **3 WaLO-OIIO 

BCU U (XLtXR»TL»YR»YLP',YRP) = (3.*(YL-vR)-(YLP'*2.«Y9P)*(XL-Xfi) ) / < XL-XRW aLO 7 1 20 

-W a LO 01 3 0 


YCU? (XP»XW« YR f YRPt AA'.QB) sAA* (xP-xR) * # 3*eB* (XP-XR) *«2-» YRP« <XP-XR) 
1 ♦ YR ■ 

Y P CURLX P . XP_»_YRR-« aA-»-8b4-s3-» * AA#4XP-XdL« «2*2,oPB« 4-X R- XR-M-Y« P — 

1000 FORMAT(20AA) 

1001 FORMAT (3P1 *4,3110) 

4002- FORM AT (AX *1 A HGEOmETR v—I-N P U-T-MX-* ZHtc* f-7-. 3tS X, 5 HHB A«a« F7, 4 , 5 X 

1 • 6HRn8ARs » F ( ,4 ) 

1003 F0R M AT (4X*BHC0ViL L IP /5X ,4HXBAR . 6X ♦ 4MY0AR • NX . 3HYPR > 

M)0A-F0RM4T («»X» I JMEXTE-RN-Ai- C4»taU>'5X^4Hx8Afl»6x,-4HY8AR,6X-»3HYPR 

1 ,10. A *4HLFUK ) 

1005 F0RHrtT(4X»lJHINTERNAi COWL/5X,4HxBaR»6X»4HYRAR,6X *3HYPR 

1 , loX t 4HLFUM — ■ : — 

DATA LAPP/ / 

if (Camp.ne.u) go to Too 

real) <5, loan - : : 

WRITE ( 6 * 1 0 r . 0 ) A 
PI02 = PII/2. 

1 — _ — TP-I-0gs3, «P 1 02 — ^ — 

READ (5* 1 00 ) ) ELL*HbADfRNBAR*JN0S8».j0UTB*JINB 
WRITE (6* I 0'2) ELL»H8 Ao*RM0aR 

RN=Rts.BAR/ELL_ : 1 : : : — 

R0=3.*PN 

RD=No-RN 

^HHaHBAR/ELI ; __ 

XFai ,*«U 

DUM=RMFLO*EM*( (1,»G0«Em**2)/gE)»a(-GC/2.) 

tt=Hh*dum : _ 

IF(LA.EU.l) TT=HH*SOnT »nUM) 


W aLO 7 140 
WaLO* 150 
—-WaLOO-160 
WAL07170 
WAL00180 
t-WaL0 7 190 
, WALC7200 
WAL00210 
-WAt-07220 
WAL07230 
WAL07240 
— WaL 0')250 
WAL00260 
WAU07270 
— w aLO ‘1280 
WAL07290 
WAL07300 
-HAL0O310 
WAL00320 
WaLO 7330 
— W aLO-7 340 
W aLO 0 350 
WaLO n 360 
-4i AL0O370 
WaL0*380 

W AL07390 
-WAL09400 
WAL07410 


EP=*000'.il 

_C ..COWL LIP- DATA — i 

JNOSajNOSb . 

IF < JNOS.EO.O) GO TC p5 

-WRI Te (6« 10:3) 

DO Io J = 1 * jN05 

REAL) <5, 1 001 ) XBAR«YHaR, YPBAR 


WAL00420 
— W AL0O430 
■ W aLO 7 44 0 
WALO 0450 
~WAL-0f>460 
WAL07470 
k«l nnitin 
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JXRI TEL16»1i) Q 1 ).XB ARtYR 4 .fl, YPBAR_ 
XN0?3 XHAR/EUl*R0 
YNOSsYHAR/ELl 

-IFXJ.EQ .1 L. p O .TO. _5 

IF (U.EO.JNCb) GO 


THElj»= p H' 
J30— TC-JL 


T C x 

a TAN { <YNOc-HH)/<RO-XNOSl ) 


toAL0n49f» 

HAL00500 

toAt0o510 

-toAL00520 

toAL0n530 

toAl0n540 

KAL00550 


— S — THE.(j ).=P-IO£ : 

GO TC 7 

6 THE < J) =TPIC2 

— 7 1 F < The ( U > • G T . p 1 1 -Er . *N(> , T HE W}-, LT . O I I+EP )— THE T^y -Pl I 

BNO^(J) =SJ p T( (XNOS-R')*«2*<YN0S-HH)*«2) 

0PNO?(J)=O. 

if <Tre i j> .Kt.Pi i >bpm/»S(j-)-3(yprar *4 mh - xmos >-/ (SIni-the ( j> > < 

ithe ( j) ) - cos ( the ( jn «opbar) i -bnos i j>*cos (The ( j > ) > /sinithe < j ) > 

lo CONTINUE 
JNC>1=0N0S-J. 


00 2<) J=1 ♦ JN0S1 

AN 1 ( j ) = ACUF ( TmE ( J) * TuE ( J* 1 ) * BNOS ( J) ♦ BNOS ( J* 1 > « BPN'Os ( J ) ,8pN0S<J*1 ) 

2n__ AN2 Cjn=BCOa-(TME ( J>-*T«E4J*-VT-+BN0S(^v-rBN0S4-^H-»fiPN0SK7) tBPNoSEJ*! ) 

C EXTERNAL CC*L data 
25 JOUTa JOUTd* 1 
JtOU I ( 44 =0. 


YOut(l)-HH 
YPCOT < 1 ) = 0. 
-LFU&C-l U =3- 


VH R I T e ( to * 1 0 ' «► ) 

DO 3fl J=2«jOUT 

HEARTS i.IC v‘.l J._.XBARi.Yai>.Rt.YP-OU-l(.JI «_lXuNO (-J-I- 

>RI.Te (6. 1 Of l ) *BAR, Y9<\R, YPOUT (J) » lFhNO < J) 
XOui < J> =XfctiH/ELL*RO 

3r. _youT.( ji =Yd aR/ell. __ - 


YEsYCUT ( JOL I ) 

JOuTisJOUT- i 

— o o -4 6:~J = 1 1 -Out i __ 

IF (LFUNO(J) • EO • 1 ) GO To 35 

AOl ( J)=ACUG (XOUT ( J) » vOtJT ( J* l) • Y 0 • 1 T ’< J I *YOUT ( J*1 ) t YPoUT ( J ) ♦ YPOUT ( J + 

Jll 


-WAL04S60 
WAL00570 
toALOoSBO 
-H*At-0f)590 
walo^goo 
toAl0O6l0 
*»aL0-1620 
1»aL00630 
WAL0-1640 
— I* A L 0 4 6 5 0 
l»AL0^660 
) WAL00670 
) WAL0T680 
V»AL00690 
(XAL00700 
-toAL0O7l0 
toAL0n7?0 
toAL0n730 
-WAL0O740 
to ALO o750 
toAL0n760 
J»AL0'i770 
to aLO 1780 
to AL On 790 
-toALOnflOO 
toALOnflIO 
toAL0082C 
— to ALOnft30 
toAL0n840 
1 to a L 0 n 8 5 0 
-toALO i860 


AO? (.J) sBCUB (XOUT (J) *vOuT(j«l) »YOljT f J) » YOUT ( J* \ ) » Y p 0UT ( J) ♦ YPOUT (J*l to aL 0 0870 

to aLO n880 

to AL0O890 


1) ) 

_G0_.Tc_ 4i! 


35 AO I ( J) si; , 
AO? (j)=;:. 
.4J. 1 -CONTINUE 


INTtfiNAL CCtoL 
IF ( JINHiEO.O) 
_JlN=jlNd>3 

JlNbl=Jl.Nd*i 
J I NB? = J I Nd *2 

_XJ Nil.) = iL. 

YlN ‘ 1 ) =HH 
YP:IN(I )=3. 

LFliN i < 1 ) =3 

WRITE (ft * 100t>) 

00 50 J=2*jlNRl 


oata 
GO TC 
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WALOn9O0 
to ACO n9 1 0 
-to AL 00920 
to aLO n930 
to aLO n940 
-to ALOn950 
toAL00960 
WAL00970 
-to ALO 0980 
toAL00990 
WALO 1000 
-to aL 0 1 0 1 0 
WAL01020 
to AUO 1 630 



— KfcAU (*• A (AU.J.J XP4R » TH *LFVvI-UJ) W AL01 0*0 

V»RITe ( 6«100i) xeaR.YB,R,YPIN(J) ,LFUmI { J) WflL01o50 

XIN <J> =*eA-/ELL*RD W aLO l 060 

Sa-JHALtjJ =Y9A"/ELL *ALOK>70 

XCbX I N ( J I N‘: 1 ) *AL01o80 

XlN(jlN«2) = (XF*XC) /2. WAL01r>90 

7 iN-(jlNbZ).= II. MaLOI 1 OO 


YPlN <JINB2)»3. 


LFllNt ( JlNd? > =3 
XlN(jIN) =XF 
_Y_lN-(jlNl sHR 


YPI'.M JIM> = 

LFU^I ( JIN) sJ 
0 0 _ g 0_J s 1. U 1 N« 2- 


IF(LFUNKJ) •E0.1)6C rO 55 

A 1 1 ( J ) =ACU9 ( X I N C J ) *XtN(J* 1) * Y I N ( J ) « Y I N ( J ♦ 1 ) t Y P I N (J) »YPIN(J*1) ) 
JU 2-1 J l = b cua ( X l N 1^1* X r iV-t J *1-1 -t-Y-lNXJ l , YP-IK-U ) * YP I N ( J ♦ 1 1 > 


-WALOl-1 10 
WAL0U20 
WALO 1130 
-*Al01 140 
WALOU50 
WAL01 160 
WaLOI 170 

WAL0U80 
W aLO 1 1 90 
W aLq l 200 


GO TO 6i' 

55 All lj)=0. 

AJ2 ( 

6" CONTINUE 
65 lAM^s 1 

10 0— I E-(..Ce. .EU.21- -GO— tO— 4 0 * 

IF (Lrfg.E^.«) GO TC ->00 

IF(LREGiGT. l.ANO.LKEri.LT.6) GO TO ?00 

C EXT Eg n a L -CC W L_rATA — 

IF {X.lT.O. )t>0 TO 120 
IF(A.c,T.xFlGo TO 130 

00-l-l.ft~Js2.« wlOuT- — 


L = J 

IF (X.lT.XOLT (J) )GO Tn 115 

Ua_COnIINUE- 

115 Y = yClH(X,XCUT(L) * YCUt ( L ) • YPOUT ( L ) »«O l <L-1> .A02(L-1) ) 
YPa^PCUb ( X • *01 JT (L 1 * Y pOu T (L 1 *A01 (L"l 1 *A02(L-1 ) ) 

: RETURN — 

1 20 YshH 
YP = ., 

— RETURN : 

130 Y=yt 

yp=o. 

: — RETURN : 

C INTERNAL CCWL GEOMEToY 

200 IF<X.LT.O«-UR.X.GT.Xr) GO TO 22 0 

00— 24-0— J=2 > J I N 

L»J 

IF<X»LT.XlMj))GO TO 215 

21-0— CON LlNUE- -^ : — 

215 YsYClR (X,XIN( L ) ,YlMi 1 .YPIMLI tAIl (L-l) *AT2(L-1> ) 
YPaYPCUW(X,Xl N (L) ,vp^N(L» tAIl (L-J ) .Al2(L-1 1 1 

RETURN : — 

220 Y=Hb 


* P«. • 

-RETUrN- 


C COvL LIP GEOMETRY 
300 DO f 1 0 J=2 • JnoS 


WaLO 1 2 1 0 
WAL01220 
— WAL-0 1 ?30 
WAL01240 
WALO 1250 
—WALOI 260 
WAL01270 
WAL01280 
-WALOI 290 
h ALO 1 300 
WAL01310 
-*AL01320 
WAL01330 
maLO 13*0 
WaLOI 350 
WaLO 1360 
WAL01370 
-WAL01380 
WAL01390 
WaLOUOO 
- walohio 

WAL01420 
WALO 1430 
-ViA L 0 1440 
WAL01450 
WAL01460 
-WAL01470 
WaLO 1480 
WAL01490 
-WaLO 1500 
WaLO 1510 
WaLO 1520 
—WaLO 1 530 
W ALO 1540 
WaLO 1 550 
-WaLO 1560 
WAL01570 
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La j. 

IF(*.LT,THP ( J) ) GO TO 315 ' 

31 0 CONTINUE 

31.5_¥ = y£uB (X « Ttt’iU) .0NCS f LL.BP.NOSiL»-*AM.l-CL-l-L.AN2<L^l U 

YP = YPCUb(X, THE <L) »6 PmOS(U »ANl (L-l) »AN2 tL-1 ) ) 

RETURN 

~C -CENT E RH 00 Y-CeokE-TR Y— — 


-a ALO 1 590 
MAL01600 
a ALO 1610 
*i«ALO 1620 
to ALO 1 630 
a AlO l 640 
a ALO l 650 


A0C-Y=0 — — 
YPaO 
RETURN 
END- 


■a ALO 1660 
aAL0l670 

a alo i6flo 

arAL0l690 


SUrRcUT InE -bOi)NOfZ«PV- 

COMpCN/bLKf!''PRATlO*RnIST«PRAD*KPLUuE* JJ 
C0MMCN/BLK6/ 7K ( 200 1 .PK < 200 ) »KEND 

Vp l — FORMAT < 2F lx .4 )- 

102 FORM A T ( 1 X ♦ I ^ * 3E 16 • A ) 

30 DO 40 K=2»«END 

! F <-2 K ( 2>— GO— Ta-50 

40 continue 

50 L1»K-1 ' 

: L2 S K - •• ' — : 

\ EP*(Z-ZK <Ll) )/ <ZK (La)-ZK(Ll ) ) 

PrES=pk(L 1 ) *EP4 ( PK (i 2) -PK (LI ) ) 

P-sifLOG (PR&VPRAT IO) 

Rf Turn 
enO 


■BoUO rt O 1 o— 

BOU0TO20 

BoUO'>o30 

BoU0l04O — 

BQU00050 

BoUOno60 

■BoUO TOZ©— 

BO^0T080 

80U00090 

■bouotioo- 

BOUOOUO 
BOU01120 
SOUO-»130" 
BOUO 1140 
BOUO'JISU 


